You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

ztbsv_L.c 4.7 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143
  1. /*********************************************************************/
  2. /* Copyright 2009, 2010 The University of Texas at Austin. */
  3. /* All rights reserved. */
  4. /* */
  5. /* Redistribution and use in source and binary forms, with or */
  6. /* without modification, are permitted provided that the following */
  7. /* conditions are met: */
  8. /* */
  9. /* 1. Redistributions of source code must retain the above */
  10. /* copyright notice, this list of conditions and the following */
  11. /* disclaimer. */
  12. /* */
  13. /* 2. Redistributions in binary form must reproduce the above */
  14. /* copyright notice, this list of conditions and the following */
  15. /* disclaimer in the documentation and/or other materials */
  16. /* provided with the distribution. */
  17. /* */
  18. /* THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF TEXAS AT */
  19. /* AUSTIN ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, */
  20. /* INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF */
  21. /* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
  22. /* DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF TEXAS AT */
  23. /* AUSTIN OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
  24. /* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES */
  25. /* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE */
  26. /* GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR */
  27. /* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF */
  28. /* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT */
  29. /* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT */
  30. /* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
  31. /* POSSIBILITY OF SUCH DAMAGE. */
  32. /* */
  33. /* The views and conclusions contained in the software and */
  34. /* documentation are those of the authors and should not be */
  35. /* interpreted as representing official policies, either expressed */
  36. /* or implied, of The University of Texas at Austin. */
  37. /*********************************************************************/
  38. #include <stdio.h>
  39. #include <ctype.h>
  40. #include "common.h"
  41. // const static FLOAT dp1 = 1.;
  42. int CNAME(BLASLONG n, BLASLONG k, FLOAT *a, BLASLONG lda, FLOAT *b, BLASLONG incb, void *buffer){
  43. BLASLONG i;
  44. FLOAT *B = b;
  45. BLASLONG length;
  46. #if (TRANSA == 2) || (TRANSA == 4)
  47. OPENBLAS_COMPLEX_FLOAT temp;
  48. #endif
  49. #ifndef UNIT
  50. FLOAT ar, ai, br, bi, ratio, den;
  51. #endif
  52. if (incb != 1) {
  53. B = buffer;
  54. COPY_K(n, b, incb, buffer, 1);
  55. }
  56. for (i = 0; i < n; i++) {
  57. #if (TRANSA == 2) || (TRANSA == 4)
  58. length = i;
  59. if (length > k) length = k;
  60. if (length > 0) {
  61. #if TRANSA == 2
  62. temp = DOTU_K(length, a + (k - length) * COMPSIZE, 1, B + (i - length) * COMPSIZE, 1);
  63. #else
  64. temp = DOTC_K(length, a + (k - length) * COMPSIZE, 1, B + (i - length) * COMPSIZE, 1);
  65. #endif
  66. B[i * 2 + 0] -= CREAL(temp);
  67. B[i * 2 + 1] -= CIMAG(temp);
  68. }
  69. #endif
  70. #ifndef UNIT
  71. #if (TRANSA == 1) || (TRANSA == 3)
  72. ar = a[0];
  73. ai = a[1];
  74. #else
  75. ar = a[k * 2 + 0];
  76. ai = a[k * 2 + 1];
  77. #endif
  78. if (fabs(ar) >= fabs(ai)){
  79. ratio = ai / ar;
  80. den = 1./(ar * ( 1 + ratio * ratio));
  81. ar = den;
  82. #if TRANSA < 3
  83. ai = -ratio * den;
  84. #else
  85. ai = ratio * den;
  86. #endif
  87. } else {
  88. ratio = ar / ai;
  89. den = 1./(ai * ( 1 + ratio * ratio));
  90. ar = ratio * den;
  91. #if TRANSA < 3
  92. ai = -den;
  93. #else
  94. ai = den;
  95. #endif
  96. }
  97. br = B[i * 2 + 0];
  98. bi = B[i * 2 + 1];
  99. B[i * 2 + 0] = ar*br - ai*bi;
  100. B[i * 2 + 1] = ar*bi + ai*br;
  101. #endif
  102. #if (TRANSA == 1) || (TRANSA == 3)
  103. length = n - i - 1;
  104. if (length > k) length = k;
  105. if (length > 0) {
  106. #if TRANSA == 1
  107. AXPYU_K(length, 0, 0,
  108. -B[i * 2 + 0], -B[i * 2 + 1],
  109. a + COMPSIZE, 1, B + (i + 1) * COMPSIZE, 1, NULL, 0);
  110. #else
  111. AXPYC_K(length, 0, 0,
  112. -B[i * 2 + 0], -B[i * 2 + 1],
  113. a + COMPSIZE, 1, B + (i + 1) * COMPSIZE, 1, NULL, 0);
  114. #endif
  115. }
  116. #endif
  117. a += lda * COMPSIZE;
  118. }
  119. if (incb != 1) {
  120. COPY_K(n, buffer, 1, b, incb);
  121. }
  122. return 0;
  123. }