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
39 #include <stdio.h>
40 #include <ctype.h>
41 #include "common.h"
42
43 const static FLOAT dm1 = -1.;
44
CNAME(BLASLONG m,FLOAT * a,FLOAT * b,BLASLONG incb,void * buffer)45 int CNAME(BLASLONG m, FLOAT *a, FLOAT *b, BLASLONG incb, void *buffer){
46
47 BLASLONG i;
48 #if (TRANSA == 2) || (TRANSA == 4)
49 FLOAT _Complex result;
50 #endif
51 #ifndef UNIT
52 FLOAT ar, ai, br, bi, ratio, den;
53 #endif
54 FLOAT *gemvbuffer = (FLOAT *)buffer;
55 FLOAT *B = b;
56
57 if (incb != 1) {
58 B = buffer;
59 gemvbuffer = (FLOAT *)(((BLASLONG)buffer + m * sizeof(FLOAT) * 2 + 4095) & ~4095);
60 COPY_K(m, b, incb, buffer, 1);
61 }
62
63 for (i = 0; i < m; i++) {
64
65 #if (TRANSA == 2) || (TRANSA == 4)
66 if (i > 0) {
67 #if TRANSA == 2
68 result = DOTU_K(i, a, 1, B, 1);
69 #else
70 result = DOTC_K(i, a, 1, B, 1);
71 #endif
72
73 B[i * COMPSIZE + 0] -= CREAL(result);
74 B[i * COMPSIZE + 1] -= CIMAG(result);
75 }
76 #endif
77
78 #ifndef UNIT
79 #if (TRANSA == 1) || (TRANSA == 3)
80 ar = a[0];
81 ai = a[1];
82 #else
83 ar = a[i * COMPSIZE + 0];
84 ai = a[i * COMPSIZE + 1];
85 #endif
86
87 if (fabs(ar) >= fabs(ai)){
88 ratio = ai / ar;
89 den = 1./(ar * ( 1 + ratio * ratio));
90
91 ar = den;
92 #if TRANSA < 3
93 ai = -ratio * den;
94 #else
95 ai = ratio * den;
96 #endif
97 } else {
98 ratio = ar / ai;
99 den = 1./(ai * ( 1 + ratio * ratio));
100 ar = ratio * den;
101 #if TRANSA < 3
102 ai = -den;
103 #else
104 ai = den;
105 #endif
106 }
107
108 br = B[i * COMPSIZE + 0];
109 bi = B[i * COMPSIZE + 1];
110
111 B[i * COMPSIZE + 0] = ar*br - ai*bi;
112 B[i * COMPSIZE + 1] = ar*bi + ai*br;
113 #endif
114
115 #if (TRANSA == 1) || (TRANSA == 3)
116 if (i < m - 1) {
117 #if TRANSA == 1
118 AXPYU_K(m - i - 1 , 0, 0,
119 - B[i * COMPSIZE + 0], - B[i * COMPSIZE + 1],
120 a + COMPSIZE, 1, B + (i + 1) * COMPSIZE, 1, NULL, 0);
121 #else
122 AXPYC_K(m - i - 1 , 0, 0,
123 - B[i * COMPSIZE + 0], - B[i * COMPSIZE + 1],
124 a + COMPSIZE, 1, B + (i + 1) * COMPSIZE, 1, NULL, 0);
125 #endif
126 }
127 #endif
128
129 #if (TRANSA == 1) || (TRANSA == 3)
130 a += (m - i) * 2;
131 #else
132 a += (i + 1) * 2;
133 #endif
134 }
135
136 if (incb != 1) {
137 COPY_K(m, buffer, 1, b, incb);
138 }
139
140 return 0;
141 }
142
143