1 /*********************************************************************/
2 /* */
3 /* Optimized BLAS libraries */
4 /* By Kazushige Goto <kgoto@tacc.utexas.edu> */
5 /* */
6 /* Copyright (c) The University of Texas, 2009. All rights reserved. */
7 /* UNIVERSITY EXPRESSLY DISCLAIMS ANY AND ALL WARRANTIES CONCERNING */
8 /* THIS SOFTWARE AND DOCUMENTATION, INCLUDING ANY WARRANTIES OF */
9 /* MERCHANTABILITY, FITNESS FOR ANY PARTICULAR PURPOSE, */
10 /* NON-INFRINGEMENT AND WARRANTIES OF PERFORMANCE, AND ANY WARRANTY */
11 /* THAT MIGHT OTHERWISE ARISE FROM COURSE OF DEALING OR USAGE OF */
12 /* TRADE. NO WARRANTY IS EITHER EXPRESS OR IMPLIED WITH RESPECT TO */
13 /* THE USE OF THE SOFTWARE OR DOCUMENTATION. */
14 /* Under no circumstances shall University be liable for incidental, */
15 /* special, indirect, direct or consequential damages or loss of */
16 /* profits, interruption of business, or related expenses which may */
17 /* arise from use of Software or Documentation, including but not */
18 /* limited to those resulting from defects in Software and/or */
19 /* Documentation, or loss or inaccuracy of data of any kind. */
20 /*********************************************************************/
21
22 #include <stdio.h>
23 #include "common.h"
24
25 #ifdef UNIT
26 #define ZTRMV ZTRMV_NLU
27 #else
28 #define ZTRMV ZTRMV_NLN
29 #endif
30
CNAME(blas_arg_t * args,BLASLONG * range_m,BLASLONG * range_n,FLOAT * sa,FLOAT * sb,BLASLONG myid)31 blasint CNAME(blas_arg_t *args, BLASLONG *range_m, BLASLONG *range_n, FLOAT *sa, FLOAT *sb, BLASLONG myid) {
32
33 BLASLONG n, lda;
34 FLOAT *a;
35
36 FLOAT ajj_r, ajj_i;
37 #ifndef UNIT
38 FLOAT ratio, den;
39 #endif
40 BLASLONG j;
41
42 n = args -> n;
43 a = (FLOAT *)args -> a;
44 lda = args -> lda;
45
46 if (range_n) {
47 n = range_n[1] - range_n[0];
48 a += range_n[0] * (lda + 1) * COMPSIZE;
49 }
50
51 for (j = n - 1; j >= 0; j--) {
52
53 ajj_r = ONE;
54 ajj_i = ZERO;
55
56 #ifndef UNIT
57 ajj_r = *(a + (j + j * lda) * COMPSIZE + 0);
58 ajj_i = *(a + (j + j * lda) * COMPSIZE + 1);
59
60 if (fabs(ajj_r) >= fabs(ajj_i)){
61 ratio = ajj_i / ajj_r;
62 den = 1. / (ajj_r * ( 1 + ratio * ratio));
63 ajj_r = den;
64 ajj_i = -ratio * den;
65 } else {
66 ratio = ajj_r / ajj_i;
67 den = 1. /(ajj_i * ( 1 + ratio * ratio));
68 ajj_r = ratio * den;
69 ajj_i = -den;
70 }
71
72 *(a + (j + j * lda) * COMPSIZE + 0) = ajj_r;
73 *(a + (j + j * lda) * COMPSIZE + 1) = ajj_i;
74 #endif
75
76 ZTRMV (n - j - 1,
77 a + ((j + 1) + (j + 1) * lda) * COMPSIZE, lda,
78 a + ((j + 1) + j * lda) * COMPSIZE, 1,
79 sb);
80
81 SCAL_K(n - j - 1, 0, 0,
82 -ajj_r, -ajj_i,
83 a + ((j + 1) + j * lda) * COMPSIZE, 1,
84 NULL, 0, NULL, 0);
85 }
86
87 return 0;
88 }
89