1
2 /*
3 This is included by sbaij.c to generate unsigned short and regular versions of these two functions
4 */
5
6 /* We cut-and-past below from aij.h to make a "no_function" version of PetscSparseDensePlusDot().
7 * This is necessary because the USESHORT case cannot use the inlined functions that may be employed. */
8
9 #if defined(PETSC_KERNEL_USE_UNROLL_4)
10 #define PetscSparseDensePlusDot_no_function(sum,r,xv,xi,nnz) { \
11 if (nnz > 0) { \
12 PetscInt nnz2=nnz,rem=nnz&0x3; \
13 switch (rem) { \
14 case 3: sum += *xv++ *r[*xi++]; \
15 case 2: sum += *xv++ *r[*xi++]; \
16 case 1: sum += *xv++ *r[*xi++]; \
17 nnz2 -= rem;} \
18 while (nnz2 > 0) { \
19 sum += xv[0] * r[xi[0]] + xv[1] * r[xi[1]] + \
20 xv[2] * r[xi[2]] + xv[3] * r[xi[3]]; \
21 xv += 4; xi += 4; nnz2 -= 4; \
22 } \
23 xv -= nnz; xi -= nnz; \
24 } \
25 }
26
27 #elif defined(PETSC_KERNEL_USE_UNROLL_2)
28 #define PetscSparseDensePlusDot_no_function(sum,r,xv,xi,nnz) { \
29 PetscInt __i,__i1,__i2; \
30 for (__i=0; __i<nnz-1; __i+=2) {__i1 = xi[__i]; __i2=xi[__i+1]; \
31 sum += (xv[__i]*r[__i1] + xv[__i+1]*r[__i2]);} \
32 if (nnz & 0x1) sum += xv[__i] * r[xi[__i]];}
33
34 #else
35 #define PetscSparseDensePlusDot_no_function(sum,r,xv,xi,nnz) { \
36 PetscInt __i; \
37 for (__i=0; __i<nnz; __i++) sum += xv[__i] * r[xi[__i]];}
38 #endif
39
40
41 #if defined(USESHORT)
MatMult_SeqSBAIJ_1_ushort(Mat A,Vec xx,Vec zz)42 PetscErrorCode MatMult_SeqSBAIJ_1_ushort(Mat A,Vec xx,Vec zz)
43 #else
44 PetscErrorCode MatMult_SeqSBAIJ_1(Mat A,Vec xx,Vec zz)
45 #endif
46 {
47 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
48 const PetscScalar *x;
49 PetscScalar *z,x1,sum;
50 const MatScalar *v;
51 MatScalar vj;
52 PetscErrorCode ierr;
53 PetscInt mbs=a->mbs,i,j,nz;
54 const PetscInt *ai=a->i;
55 #if defined(USESHORT)
56 const unsigned short *ib=a->jshort;
57 unsigned short ibt;
58 #else
59 const PetscInt *ib=a->j;
60 PetscInt ibt;
61 #endif
62 PetscInt nonzerorow=0,jmin;
63 #if defined(PETSC_USE_COMPLEX)
64 const int aconj = A->hermitian;
65 #else
66 const int aconj = 0;
67 #endif
68
69 PetscFunctionBegin;
70 ierr = VecSet(zz,0.0);CHKERRQ(ierr);
71 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr);
72 ierr = VecGetArray(zz,&z);CHKERRQ(ierr);
73
74 v = a->a;
75 for (i=0; i<mbs; i++) {
76 nz = ai[i+1] - ai[i]; /* length of i_th row of A */
77 if (!nz) continue; /* Move to the next row if the current row is empty */
78 nonzerorow++;
79 sum = 0.0;
80 jmin = 0;
81 x1 = x[i];
82 if (ib[0] == i) {
83 sum = v[0]*x1; /* diagonal term */
84 jmin++;
85 }
86 PetscPrefetchBlock(ib+nz,nz,0,PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
87 PetscPrefetchBlock(v+nz,nz,0,PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
88 if (aconj) {
89 for (j=jmin; j<nz; j++) {
90 ibt = ib[j];
91 vj = v[j];
92 z[ibt] += PetscConj(vj) * x1; /* (strict lower triangular part of A)*x */
93 sum += vj * x[ibt]; /* (strict upper triangular part of A)*x */
94 }
95 } else {
96 for (j=jmin; j<nz; j++) {
97 ibt = ib[j];
98 vj = v[j];
99 z[ibt] += vj * x1; /* (strict lower triangular part of A)*x */
100 sum += vj * x[ibt]; /* (strict upper triangular part of A)*x */
101 }
102 }
103 z[i] += sum;
104 v += nz;
105 ib += nz;
106 }
107
108 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr);
109 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr);
110 ierr = PetscLogFlops(2.0*(2.0*a->nz - nonzerorow) - nonzerorow);CHKERRQ(ierr);
111 PetscFunctionReturn(0);
112 }
113
114 #if defined(USESHORT)
MatSOR_SeqSBAIJ_ushort(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)115 PetscErrorCode MatSOR_SeqSBAIJ_ushort(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
116 #else
117 PetscErrorCode MatSOR_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
118 #endif
119 {
120 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
121 const MatScalar *aa=a->a,*v,*v1,*aidiag;
122 PetscScalar *x,*t,sum;
123 const PetscScalar *b;
124 MatScalar tmp;
125 PetscErrorCode ierr;
126 PetscInt m =a->mbs,bs=A->rmap->bs,j;
127 const PetscInt *ai=a->i;
128 #if defined(USESHORT)
129 const unsigned short *aj=a->jshort,*vj,*vj1;
130 #else
131 const PetscInt *aj=a->j,*vj,*vj1;
132 #endif
133 PetscInt nz,nz1,i;
134
135 PetscFunctionBegin;
136 if (fshift == -1.0) fshift = 0.0; /* negative fshift indicates do not error on zero diagonal; this code never errors on zero diagonal */
137 if (flag & SOR_EISENSTAT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for Eisenstat");
138
139 its = its*lits;
140 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
141
142 if (bs > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented");
143
144 ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
145 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr);
146
147 if (!a->idiagvalid) {
148 if (!a->idiag) {
149 ierr = PetscMalloc1(m,&a->idiag);CHKERRQ(ierr);
150 }
151 for (i=0; i<a->mbs; i++) a->idiag[i] = 1.0/a->a[a->i[i]];
152 a->idiagvalid = PETSC_TRUE;
153 }
154
155 if (!a->sor_work) {
156 ierr = PetscMalloc1(m,&a->sor_work);CHKERRQ(ierr);
157 }
158 t = a->sor_work;
159
160 aidiag = a->idiag;
161
162 if (flag == SOR_APPLY_UPPER) {
163 /* apply (U + D/omega) to the vector */
164 PetscScalar d;
165 for (i=0; i<m; i++) {
166 d = fshift + aa[ai[i]];
167 nz = ai[i+1] - ai[i] - 1;
168 vj = aj + ai[i] + 1;
169 v = aa + ai[i] + 1;
170 sum = b[i]*d/omega;
171 #ifdef USESHORT
172 PetscSparseDensePlusDot_no_function(sum,b,v,vj,nz);
173 #else
174 PetscSparseDensePlusDot(sum,b,v,vj,nz);
175 #endif
176 x[i] = sum;
177 }
178 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr);
179 }
180
181 if (flag & SOR_ZERO_INITIAL_GUESS) {
182 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
183 ierr = PetscArraycpy(t,b,m);CHKERRQ(ierr);
184
185 v = aa + 1;
186 vj = aj + 1;
187 for (i=0; i<m; i++) {
188 nz = ai[i+1] - ai[i] - 1;
189 tmp = -(x[i] = omega*t[i]*aidiag[i]);
190 for (j=0; j<nz; j++) t[vj[j]] += tmp*v[j];
191 v += nz + 1;
192 vj += nz + 1;
193 }
194 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr);
195 }
196
197 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
198 int nz2;
199 if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)) {
200 #if defined(PETSC_USE_BACKWARD_LOOP)
201 v = aa + ai[m] - 1;
202 vj = aj + ai[m] - 1;
203 for (i=m-1; i>=0; i--) {
204 sum = b[i];
205 nz = ai[i+1] - ai[i] - 1;
206 {PetscInt __i;for (__i=0; __i<nz; __i++) sum -= v[-__i] * x[vj[-__i]];}
207 #else
208 v = aa + ai[m-1] + 1;
209 vj = aj + ai[m-1] + 1;
210 nz = 0;
211 for (i=m-1; i>=0; i--) {
212 sum = b[i];
213 nz2 = ai[i] - ai[PetscMax(i-1,0)] - 1; /* avoid referencing ai[-1], nonsense nz2 is okay on last iteration */
214 PETSC_Prefetch(v-nz2-1,0,PETSC_PREFETCH_HINT_NTA);
215 PETSC_Prefetch(vj-nz2-1,0,PETSC_PREFETCH_HINT_NTA);
216 PetscSparseDenseMinusDot(sum,x,v,vj,nz);
217 nz = nz2;
218 #endif
219 x[i] = omega*sum*aidiag[i];
220 v -= nz + 1;
221 vj -= nz + 1;
222 }
223 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr);
224 } else {
225 v = aa + ai[m-1] + 1;
226 vj = aj + ai[m-1] + 1;
227 nz = 0;
228 for (i=m-1; i>=0; i--) {
229 sum = t[i];
230 nz2 = ai[i] - ai[PetscMax(i-1,0)] - 1; /* avoid referencing ai[-1], nonsense nz2 is okay on last iteration */
231 PETSC_Prefetch(v-nz2-1,0,PETSC_PREFETCH_HINT_NTA);
232 PETSC_Prefetch(vj-nz2-1,0,PETSC_PREFETCH_HINT_NTA);
233 PetscSparseDenseMinusDot(sum,x,v,vj,nz);
234 x[i] = (1-omega)*x[i] + omega*sum*aidiag[i];
235 nz = nz2;
236 v -= nz + 1;
237 vj -= nz + 1;
238 }
239 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr);
240 }
241 }
242 its--;
243 }
244
245 while (its--) {
246 /*
247 forward sweep:
248 for i=0,...,m-1:
249 sum[i] = (b[i] - U(i,:)x)/d[i];
250 x[i] = (1-omega)x[i] + omega*sum[i];
251 b = b - x[i]*U^T(i,:);
252
253 */
254 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
255 ierr = PetscArraycpy(t,b,m);CHKERRQ(ierr);
256
257 for (i=0; i<m; i++) {
258 v = aa + ai[i] + 1; v1=v;
259 vj = aj + ai[i] + 1; vj1=vj;
260 nz = ai[i+1] - ai[i] - 1; nz1=nz;
261 sum = t[i];
262 while (nz1--) sum -= (*v1++)*x[*vj1++];
263 x[i] = (1-omega)*x[i] + omega*sum*aidiag[i];
264 while (nz--) t[*vj++] -= x[i]*(*v++);
265 }
266 ierr = PetscLogFlops(4.0*a->nz);CHKERRQ(ierr);
267 }
268
269 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
270 /*
271 backward sweep:
272 b = b - x[i]*U^T(i,:), i=0,...,n-2
273 for i=m-1,...,0:
274 sum[i] = (b[i] - U(i,:)x)/d[i];
275 x[i] = (1-omega)x[i] + omega*sum[i];
276 */
277 /* if there was a forward sweep done above then I thing the next two for loops are not needed */
278 ierr = PetscArraycpy(t,b,m);CHKERRQ(ierr);
279
280 for (i=0; i<m-1; i++) { /* update rhs */
281 v = aa + ai[i] + 1;
282 vj = aj + ai[i] + 1;
283 nz = ai[i+1] - ai[i] - 1;
284 while (nz--) t[*vj++] -= x[i]*(*v++);
285 }
286 ierr = PetscLogFlops(2.0*(a->nz - m));CHKERRQ(ierr);
287 for (i=m-1; i>=0; i--) {
288 v = aa + ai[i] + 1;
289 vj = aj + ai[i] + 1;
290 nz = ai[i+1] - ai[i] - 1;
291 sum = t[i];
292 while (nz--) sum -= x[*vj++]*(*v++);
293 x[i] = (1-omega)*x[i] + omega*sum*aidiag[i];
294 }
295 ierr = PetscLogFlops(2.0*(a->nz + m));CHKERRQ(ierr);
296 }
297 }
298
299 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
300 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr);
301 PetscFunctionReturn(0);
302 }
303