Actual source code: relax.h
2: /*
3: This is included by sbaij.c to generate unsigned short and regular versions of these two functions
4: */
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. */
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: }
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]];}
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
41: #if defined(USESHORT)
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
70: VecSet(zz,0.0);
71: VecGetArrayRead(xx,&x);
72: VecGetArray(zz,&z);
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: }
108: VecRestoreArrayRead(xx,&x);
109: VecRestoreArray(zz,&z);
110: PetscLogFlops(2.0*(2.0*a->nz - nonzerorow) - nonzerorow);
111: return(0);
112: }
114: #if defined(USESHORT)
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;
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");
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);
142: if (bs > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented");
144: VecGetArray(xx,&x);
145: VecGetArrayRead(bb,&b);
147: if (!a->idiagvalid) {
148: if (!a->idiag) {
149: PetscMalloc1(m,&a->idiag);
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: }
155: if (!a->sor_work) {
156: PetscMalloc1(m,&a->sor_work);
157: }
158: t = a->sor_work;
160: aidiag = a->idiag;
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: PetscLogFlops(a->nz);
179: }
181: if (flag & SOR_ZERO_INITIAL_GUESS) {
182: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
183: PetscArraycpy(t,b,m);
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: PetscLogFlops(2.0*a->nz);
195: }
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: PetscLogFlops(2.0*a->nz);
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: PetscLogFlops(2.0*a->nz);
240: }
241: }
242: its--;
243: }
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,:);
253: */
254: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
255: PetscArraycpy(t,b,m);
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: PetscLogFlops(4.0*a->nz);
267: }
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: PetscArraycpy(t,b,m);
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: PetscLogFlops(2.0*(a->nz - m));
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: PetscLogFlops(2.0*(a->nz + m));
296: }
297: }
299: VecRestoreArray(xx,&x);
300: VecRestoreArrayRead(bb,&b);
301: return(0);
302: }