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.

ztpsv_U.c 4.6 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133
  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. int CNAME(BLASLONG m, FLOAT *a, FLOAT *b, BLASLONG incb, void *buffer){
  42. BLASLONG i;
  43. #if (TRANSA == 2) || (TRANSA == 4)
  44. OPENBLAS_COMPLEX_FLOAT result;
  45. #endif
  46. #ifndef UNIT
  47. FLOAT ar, ai, br, bi, ratio, den;
  48. #endif
  49. FLOAT *B = b;
  50. if (incb != 1) {
  51. B = buffer;
  52. COPY_K(m, b, incb, buffer, 1);
  53. }
  54. a += (m + 1) * m - 2;
  55. for (i = 0; i < m; i++) {
  56. #if (TRANSA == 2) || (TRANSA == 4)
  57. if (i > 0) {
  58. #if TRANSA == 2
  59. result = DOTU_K(i, a + 2, 1, B + (m - i) * 2, 1);
  60. #else
  61. result = DOTC_K(i, a + 2, 1, B + (m - i) * 2, 1);
  62. #endif
  63. B[(m - i - 1) * 2 + 0] -= CREAL(result);
  64. B[(m - i - 1) * 2 + 1] -= CIMAG(result);
  65. }
  66. #endif
  67. #ifndef UNIT
  68. ar = a[0];
  69. ai = a[1];
  70. if (fabs(ar) >= fabs(ai)){
  71. ratio = ai / ar;
  72. den = 1./(ar * ( 1 + ratio * ratio));
  73. ar = den;
  74. #if (TRANSA == 1) || (TRANSA == 2)
  75. ai = -ratio * den;
  76. #else
  77. ai = ratio * den;
  78. #endif
  79. } else {
  80. ratio = ar / ai;
  81. den = 1./(ai * ( 1 + ratio * ratio));
  82. ar = ratio * den;
  83. #if (TRANSA == 1) || (TRANSA == 2)
  84. ai = -den;
  85. #else
  86. ai = den;
  87. #endif
  88. }
  89. br = B[(m - i - 1) * 2 + 0];
  90. bi = B[(m - i - 1) * 2 + 1];
  91. B[(m - i - 1) * 2 + 0] = ar*br - ai*bi;
  92. B[(m - i - 1) * 2 + 1] = ar*bi + ai*br;
  93. #endif
  94. #if (TRANSA == 1) || (TRANSA == 3)
  95. if (i < m - 1) {
  96. #if TRANSA == 1
  97. AXPYU_K (m - i - 1, 0, 0, - B[(m - i - 1) * 2 + 0], -B[(m - i - 1) * 2 + 1],
  98. a - (m - i - 1) * COMPSIZE, 1, B, 1, NULL, 0);
  99. #else
  100. AXPYC_K (m - i - 1, 0, 0, - B[(m - i - 1) * 2 + 0], -B[(m - i - 1) * 2 + 1],
  101. a - (m - i - 1) * COMPSIZE, 1, B, 1, NULL, 0);
  102. #endif
  103. }
  104. #endif
  105. #if (TRANSA == 1) || (TRANSA == 3)
  106. a -= (m - i) * 2;
  107. #else
  108. a -= (i + 2) * 2;
  109. #endif
  110. }
  111. if (incb != 1) {
  112. COPY_K(m, buffer, 1, b, incb);
  113. }
  114. return 0;
  115. }