1 //==============================================================================
2 //=== sfmult_atxtyn_k ==========================================================
3 //==============================================================================
4
5 // y = A'*x' where x has 2, 3, or 4 rows
6
7 // compare with sfmult_atxtyt
8
9 // sfmult_AT_XT_YN_2 y = A'*x' where x is 2-by-m, and y is n-by-2
10 // sfmult_AT_XT_YN_3 y = A'*x' where x is 3-by-m, and y is n-by-3 (ldx = 4)
11 // sfmult_AT_XT_YN_4 y = A'*x' where x is 4-by-m, and y is n-by-4
12
13 #include "sfmult.h"
14
15 //==============================================================================
16 //=== sfmult_AT_XT_YN_2 ========================================================
17 //==============================================================================
18
sfmult_AT_XT_YN_2(double * Yx,double * Yz,const Int * Ap,const Int * Ai,const double * Ax,const double * Az,Int m,Int n,const double * Xx,const double * Xz,int ac,int xc,int yc)19 void sfmult_AT_XT_YN_2 // y = A'*x' x is 2-by-m, and y is n-by-2
20 (
21 // --- outputs, not initialized on input
22 double *Yx, // n-by-2
23 double *Yz, // n-by-2 if Y is complex (TO DO)
24
25 // --- inputs, not modified
26 const Int *Ap, // size n+1 column pointers
27 const Int *Ai, // size nz = Ap[n] row indices
28 const double *Ax, // size nz values
29 const double *Az, // size nz imaginary values if A is complex (TO DO)
30 Int m, // A is m-by-n
31 Int n,
32 const double *Xx, // 2-by-m
33 const double *Xz, // 2-by-m if X complex (TO DO)
34 int ac, // true: use conj(A), otherwise use A (TO DO)
35 int xc, // true: use conj(X), otherwise use X (TO DO)
36 int yc // true: compute conj(Y), otherwise compute Y (TO DO)
37 )
38 {
39 double y [2], a [4] ;
40 Int p, pend, j, i0, i1, i2 ;
41
42 p = 0 ;
43 for (j = 0 ; j < n ; j++)
44 {
45 pend = Ap [j+1] ;
46 y [0] = 0 ;
47 y [1] = 0 ;
48 switch ((pend - p) % 3)
49 {
50 case 2:
51 i0 = Ai [p] ;
52 a [0] = Ax [p] ;
53 y [0] += a [0] * Xx [2*i0 ] ;
54 y [1] += a [0] * Xx [2*i0+1] ;
55 p++ ;
56 case 1:
57 i0 = Ai [p] ;
58 a [0] = Ax [p] ;
59 y [0] += a [0] * Xx [2*i0 ] ;
60 y [1] += a [0] * Xx [2*i0+1] ;
61 p++ ;
62 case 0: ;
63 }
64 for ( ; p < pend ; p += 3)
65 {
66 i0 = Ai [p ] ;
67 i1 = Ai [p+1] ;
68 i2 = Ai [p+2] ;
69 a [0] = Ax [p ] ;
70 a [1] = Ax [p+1] ;
71 a [2] = Ax [p+2] ;
72 y [0] += a [0] * Xx [2*i0 ] ;
73 y [1] += a [0] * Xx [2*i0+1] ;
74 y [0] += a [1] * Xx [2*i1 ] ;
75 y [1] += a [1] * Xx [2*i1+1] ;
76 y [0] += a [2] * Xx [2*i2 ] ;
77 y [1] += a [2] * Xx [2*i2+1] ;
78 }
79 Yx [j ] = y [0] ;
80 Yx [j+n] = y [1] ;
81 }
82 }
83
84
85 //==============================================================================
86 //=== sfmult_AT_XT_YN_3 ========================================================
87 //==============================================================================
88
sfmult_AT_XT_YN_3(double * Yx,double * Yz,const int * Ap,const int * Ai,const double * Ax,const double * Az,int m,int n,const double * Xx,const double * Xz,int ac,int xc,int yc)89 void sfmult_AT_XT_YN_3 // y = A'*x' x is 3-by-m, and y is n-by-3 (ldx = 4)
90 (
91 // --- outputs, not initialized on input
92 double *Yx, // n-by-3
93 double *Yz, // n-by-3 if Y is complex (TO DO)
94
95 // --- inputs, not modified
96 const int *Ap, // size n+1 column pointers
97 const int *Ai, // size nz = Ap[n] row indices
98 const double *Ax, // size nz values
99 const double *Az, // size nz imaginary values if A is complex (TO DO)
100 int m, // A is m-by-n
101 int n,
102 const double *Xx, // 3-by-m
103 const double *Xz, // 3-by-m if X complex (TO DO)
104 int ac, // true: use conj(A), otherwise use A (TO DO)
105 int xc, // true: use conj(X), otherwise use X (TO DO)
106 int yc // true: compute conj(Y), otherwise compute Y (TO DO)
107 )
108 {
109 double y [4], a [2] ;
110 int p, pend, j, i0, i1 ;
111
112 p = 0 ;
113 for (j = 0 ; j < n ; j++)
114 {
115 pend = Ap [j+1] ;
116 y [0] = 0 ;
117 y [1] = 0 ;
118 y [2] = 0 ;
119 if ((pend - p) % 2)
120 {
121 i0 = Ai [p] ;
122 a [0] = Ax [p] ;
123 y [0] += a [0] * Xx [4*i0 ] ;
124 y [1] += a [0] * Xx [4*i0+1] ;
125 y [2] += a [0] * Xx [4*i0+2] ;
126 p++ ;
127 }
128 for ( ; p < pend ; p += 2)
129 {
130 i0 = Ai [p ] ;
131 i1 = Ai [p+1] ;
132 a [0] = Ax [p ] ;
133 a [1] = Ax [p+1] ;
134 y [0] += a [0] * Xx [4*i0 ] ;
135 y [1] += a [0] * Xx [4*i0+1] ;
136 y [2] += a [0] * Xx [4*i0+2] ;
137 y [0] += a [1] * Xx [4*i1 ] ;
138 y [1] += a [1] * Xx [4*i1+1] ;
139 y [2] += a [1] * Xx [4*i1+2] ;
140 }
141 Yx [j ] = y [0] ;
142 Yx [j+ n] = y [1] ;
143 Yx [j+2*n] = y [2] ;
144 }
145 }
146
147
148 //==============================================================================
149 //=== sfmult_AT_XT_YN_4 ========================================================
150 //==============================================================================
151
sfmult_AT_XT_YN_4(double * Yx,double * Yz,const int * Ap,const int * Ai,const double * Ax,const double * Az,int m,int n,const double * Xx,const double * Xz,int ac,int xc,int yc)152 void sfmult_AT_XT_YN_4 // y = A'*x' x is 4-by-m, and y is n-by-4
153 (
154 // --- outputs, not initialized on input
155 double *Yx, // n-by-4
156 double *Yz, // n-by-4 if Y is complex (TO DO)
157
158 // --- inputs, not modified
159 const int *Ap, // size n+1 column pointers
160 const int *Ai, // size nz = Ap[n] row indices
161 const double *Ax, // size nz values
162 const double *Az, // size nz imaginary values if A is complex (TO DO)
163 int m, // A is m-by-n
164 int n,
165 const double *Xx, // 4-by-m
166 const double *Xz, // 4-by-m if X complex (TO DO)
167 int ac, // true: use conj(A), otherwise use A (TO DO)
168 int xc, // true: use conj(X), otherwise use X (TO DO)
169 int yc // true: compute conj(Y), otherwise compute Y (TO DO)
170 )
171 {
172 double y [4], a ;
173 int p, pend, j, i ;
174
175 p = 0 ;
176 for (j = 0 ; j < n ; j++)
177 {
178 pend = Ap [j+1] ;
179 y [0] = 0 ;
180 y [1] = 0 ;
181 y [2] = 0 ;
182 y [3] = 0 ;
183 for ( ; p < pend ; p++)
184 {
185 i = Ai [p] ;
186 a = Ax [p] ;
187 y [0] += a * Xx [4*i ] ;
188 y [1] += a * Xx [4*i+1] ;
189 y [2] += a * Xx [4*i+2] ;
190 y [3] += a * Xx [4*i+3] ;
191 }
192 Yx [j ] = y [0] ;
193 Yx [j+ n] = y [1] ;
194 Yx [j+2*n] = y [2] ;
195 Yx [j+3*n] = y [3] ;
196 }
197 }
198