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:     { \
 12:       if (nnz > 0) { \
 13:         PetscInt nnz2 = nnz, rem = nnz & 0x3; \
 14:         switch (rem) { \
 15:         case 3: \
 16:           sum += *xv++ * r[*xi++]; \
 17:         case 2: \
 18:           sum += *xv++ * r[*xi++]; \
 19:         case 1: \
 20:           sum += *xv++ * r[*xi++]; \
 21:           nnz2 -= rem; \
 22:         } \
 23:         while (nnz2 > 0) { \
 24:           sum += xv[0] * r[xi[0]] + xv[1] * r[xi[1]] + xv[2] * r[xi[2]] + xv[3] * r[xi[3]]; \
 25:           xv += 4; \
 26:           xi += 4; \
 27:           nnz2 -= 4; \
 28:         } \
 29:         xv -= nnz; \
 30:         xi -= nnz; \
 31:       } \
 32:     }

 34: #elif defined(PETSC_KERNEL_USE_UNROLL_2)
 35:   #define PetscSparseDensePlusDot_no_function(sum, r, xv, xi, nnz) \
 36:     { \
 37:       PetscInt __i, __i1, __i2; \
 38:       for (__i = 0; __i < nnz - 1; __i += 2) { \
 39:         __i1 = xi[__i]; \
 40:         __i2 = xi[__i + 1]; \
 41:         sum += (xv[__i] * r[__i1] + xv[__i + 1] * r[__i2]); \
 42:       } \
 43:       if (nnz & 0x1) sum += xv[__i] * r[xi[__i]]; \
 44:     }

 46: #else
 47:   #define PetscSparseDensePlusDot_no_function(sum, r, xv, xi, nnz) \
 48:     { \
 49:       PetscInt __i; \
 50:       for (__i = 0; __i < nnz; __i++) sum += xv[__i] * r[xi[__i]]; \
 51:     }
 52: #endif

 54: #if defined(USESHORT)
 55: PetscErrorCode MatMult_SeqSBAIJ_1_ushort(Mat A, Vec xx, Vec zz)
 56: #else
 57: PetscErrorCode MatMult_SeqSBAIJ_1(Mat A, Vec xx, Vec zz)
 58: #endif
 59: {
 60:   Mat_SeqSBAIJ      *a = (Mat_SeqSBAIJ *)A->data;
 61:   const PetscScalar *x;
 62:   PetscScalar       *z, x1, sum;
 63:   const MatScalar   *v;
 64:   MatScalar          vj;
 65:   PetscInt           mbs = a->mbs, i, j, nz;
 66:   const PetscInt    *ai  = a->i;
 67: #if defined(USESHORT)
 68:   const unsigned short *ib = a->jshort;
 69:   unsigned short        ibt;
 70: #else
 71:   const PetscInt *ib = a->j;
 72:   PetscInt        ibt;
 73: #endif
 74:   PetscInt nonzerorow = 0, jmin;
 75: #if defined(PETSC_USE_COMPLEX)
 76:   const int aconj = A->hermitian == PETSC_BOOL3_TRUE;
 77: #else
 78:   const int       aconj = 0;
 79: #endif

 81:   PetscFunctionBegin;
 82:   PetscCall(VecSet(zz, 0.0));
 83:   PetscCall(VecGetArrayRead(xx, &x));
 84:   PetscCall(VecGetArray(zz, &z));

 86:   v = a->a;
 87:   for (i = 0; i < mbs; i++) {
 88:     nz = ai[i + 1] - ai[i]; /* length of i_th row of A */
 89:     if (!nz) continue;      /* Move to the next row if the current row is empty */
 90:     nonzerorow++;
 91:     sum  = 0.0;
 92:     jmin = 0;
 93:     x1   = x[i];
 94:     if (ib[0] == i) {
 95:       sum = v[0] * x1; /* diagonal term */
 96:       jmin++;
 97:     }
 98:     PetscPrefetchBlock(ib + nz, nz, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
 99:     PetscPrefetchBlock(v + nz, nz, 0, PETSC_PREFETCH_HINT_NTA);  /* Entries for the next row */
100:     if (aconj) {
101:       for (j = jmin; j < nz; j++) {
102:         ibt = ib[j];
103:         vj  = v[j];
104:         z[ibt] += PetscConj(vj) * x1; /* (strict lower triangular part of A)*x  */
105:         sum += vj * x[ibt];           /* (strict upper triangular part of A)*x  */
106:       }
107:     } else {
108:       for (j = jmin; j < nz; j++) {
109:         ibt = ib[j];
110:         vj  = v[j];
111:         z[ibt] += vj * x1;  /* (strict lower triangular part of A)*x  */
112:         sum += vj * x[ibt]; /* (strict upper triangular part of A)*x  */
113:       }
114:     }
115:     z[i] += sum;
116:     v += nz;
117:     ib += nz;
118:   }

120:   PetscCall(VecRestoreArrayRead(xx, &x));
121:   PetscCall(VecRestoreArray(zz, &z));
122:   PetscCall(PetscLogFlops(2.0 * (2.0 * a->nz - nonzerorow) - nonzerorow));
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: #if defined(USESHORT)
127: PetscErrorCode MatSOR_SeqSBAIJ_ushort(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
128: #else
129: PetscErrorCode MatSOR_SeqSBAIJ(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
130: #endif
131: {
132:   Mat_SeqSBAIJ      *a  = (Mat_SeqSBAIJ *)A->data;
133:   const MatScalar   *aa = a->a, *v, *v1, *aidiag;
134:   PetscScalar       *x, *t, sum;
135:   const PetscScalar *b;
136:   MatScalar          tmp;
137:   PetscInt           m = a->mbs, bs = A->rmap->bs, j;
138:   const PetscInt    *ai = a->i;
139: #if defined(USESHORT)
140:   const unsigned short *aj = a->jshort, *vj, *vj1;
141: #else
142:   const PetscInt *aj = a->j, *vj, *vj1;
143: #endif
144:   PetscInt nz, nz1, i;

146:   PetscFunctionBegin;
147:   if (fshift == -1.0) fshift = 0.0; /* negative fshift indicates do not error on zero diagonal; this code never errors on zero diagonal */
148:   PetscCheck(!(flag & SOR_EISENSTAT), PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for Eisenstat");

150:   its = its * lits;
151:   PetscCheck(its > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Relaxation requires global its %" PetscInt_FMT " and local its %" PetscInt_FMT " both positive", its, lits);

153:   PetscCheck(bs <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "SSOR for block size > 1 is not yet implemented");

155:   PetscCall(VecGetArray(xx, &x));
156:   PetscCall(VecGetArrayRead(bb, &b));

158:   if (!a->idiagvalid) {
159:     if (!a->idiag) PetscCall(PetscMalloc1(m, &a->idiag));
160:     for (i = 0; i < a->mbs; i++) a->idiag[i] = 1.0 / a->a[a->i[i]];
161:     a->idiagvalid = PETSC_TRUE;
162:   }

164:   if (!a->sor_work) PetscCall(PetscMalloc1(m, &a->sor_work));
165:   t = a->sor_work;

167:   aidiag = a->idiag;

169:   if (flag == SOR_APPLY_UPPER) {
170:     /* apply (U + D/omega) to the vector */
171:     PetscScalar d;
172:     for (i = 0; i < m; i++) {
173:       d   = fshift + aa[ai[i]];
174:       nz  = ai[i + 1] - ai[i] - 1;
175:       vj  = aj + ai[i] + 1;
176:       v   = aa + ai[i] + 1;
177:       sum = b[i] * d / omega;
178: #ifdef USESHORT
179:       PetscSparseDensePlusDot_no_function(sum, b, v, vj, nz);
180: #else
181:       PetscSparseDensePlusDot(sum, b, v, vj, nz);
182: #endif
183:       x[i] = sum;
184:     }
185:     PetscCall(PetscLogFlops(a->nz));
186:   }

188:   if (flag & SOR_ZERO_INITIAL_GUESS) {
189:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
190:       PetscCall(PetscArraycpy(t, b, m));

192:       v  = aa + 1;
193:       vj = aj + 1;
194:       for (i = 0; i < m; i++) {
195:         nz  = ai[i + 1] - ai[i] - 1;
196:         tmp = -(x[i] = omega * t[i] * aidiag[i]);
197:         for (j = 0; j < nz; j++) t[vj[j]] += tmp * v[j];
198:         v += nz + 1;
199:         vj += nz + 1;
200:       }
201:       PetscCall(PetscLogFlops(2.0 * a->nz));
202:     }

204:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
205:       int nz2;
206:       if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)) {
207: #if defined(PETSC_USE_BACKWARD_LOOP)
208:         v  = aa + ai[m] - 1;
209:         vj = aj + ai[m] - 1;
210:         for (i = m - 1; i >= 0; i--) {
211:           sum = b[i];
212:           nz  = ai[i + 1] - ai[i] - 1;
213:           {
214:             PetscInt __i;
215:             for (__i = 0; __i < nz; __i++) sum -= v[-__i] * x[vj[-__i]];
216:           }
217: #else
218:         v  = aa + ai[m - 1] + 1;
219:         vj = aj + ai[m - 1] + 1;
220:         nz = 0;
221:         for (i = m - 1; i >= 0; i--) {
222:           sum = b[i];
223:           nz2 = ai[i] - ai[PetscMax(i - 1, 0)] - 1; /* avoid referencing ai[-1], nonsense nz2 is okay on last iteration */
224:           PETSC_Prefetch(v - nz2 - 1, 0, PETSC_PREFETCH_HINT_NTA);
225:           PETSC_Prefetch(vj - nz2 - 1, 0, PETSC_PREFETCH_HINT_NTA);
226:           PetscSparseDenseMinusDot(sum, x, v, vj, nz);
227:           nz = nz2;
228: #endif
229:           x[i] = omega * sum * aidiag[i];
230:           v -= nz + 1;
231:           vj -= nz + 1;
232:         }
233:         PetscCall(PetscLogFlops(2.0 * a->nz));
234:       } else {
235:         v  = aa + ai[m - 1] + 1;
236:         vj = aj + ai[m - 1] + 1;
237:         nz = 0;
238:         for (i = m - 1; i >= 0; i--) {
239:           sum = t[i];
240:           nz2 = ai[i] - ai[PetscMax(i - 1, 0)] - 1; /* avoid referencing ai[-1], nonsense nz2 is okay on last iteration */
241:           PETSC_Prefetch(v - nz2 - 1, 0, PETSC_PREFETCH_HINT_NTA);
242:           PETSC_Prefetch(vj - nz2 - 1, 0, PETSC_PREFETCH_HINT_NTA);
243:           PetscSparseDenseMinusDot(sum, x, v, vj, nz);
244:           x[i] = (1 - omega) * x[i] + omega * sum * aidiag[i];
245:           nz   = nz2;
246:           v -= nz + 1;
247:           vj -= nz + 1;
248:         }
249:         PetscCall(PetscLogFlops(2.0 * a->nz));
250:       }
251:     }
252:     its--;
253:   }

255:   while (its--) {
256:     /*
257:        forward sweep:
258:        for i=0,...,m-1:
259:          sum[i] = (b[i] - U(i,:)x)/d[i];
260:          x[i]   = (1-omega)x[i] + omega*sum[i];
261:          b      = b - x[i]*U^T(i,:);

263:     */
264:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
265:       PetscCall(PetscArraycpy(t, b, m));

267:       for (i = 0; i < m; i++) {
268:         v   = aa + ai[i] + 1;
269:         v1  = v;
270:         vj  = aj + ai[i] + 1;
271:         vj1 = vj;
272:         nz  = ai[i + 1] - ai[i] - 1;
273:         nz1 = nz;
274:         sum = t[i];
275:         while (nz1--) sum -= (*v1++) * x[*vj1++];
276:         x[i] = (1 - omega) * x[i] + omega * sum * aidiag[i];
277:         while (nz--) t[*vj++] -= x[i] * (*v++);
278:       }
279:       PetscCall(PetscLogFlops(4.0 * a->nz));
280:     }

282:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
283:       /*
284:        backward sweep:
285:        b = b - x[i]*U^T(i,:), i=0,...,n-2
286:        for i=m-1,...,0:
287:          sum[i] = (b[i] - U(i,:)x)/d[i];
288:          x[i]   = (1-omega)x[i] + omega*sum[i];
289:       */
290:       /* if there was a forward sweep done above then I thing the next two for loops are not needed */
291:       PetscCall(PetscArraycpy(t, b, m));

293:       for (i = 0; i < m - 1; i++) { /* update rhs */
294:         v  = aa + ai[i] + 1;
295:         vj = aj + ai[i] + 1;
296:         nz = ai[i + 1] - ai[i] - 1;
297:         while (nz--) t[*vj++] -= x[i] * (*v++);
298:       }
299:       PetscCall(PetscLogFlops(2.0 * (a->nz - m)));
300:       for (i = m - 1; i >= 0; i--) {
301:         v   = aa + ai[i] + 1;
302:         vj  = aj + ai[i] + 1;
303:         nz  = ai[i + 1] - ai[i] - 1;
304:         sum = t[i];
305:         while (nz--) sum -= x[*vj++] * (*v++);
306:         x[i] = (1 - omega) * x[i] + omega * sum * aidiag[i];
307:       }
308:       PetscCall(PetscLogFlops(2.0 * (a->nz + m)));
309:     }
310:   }

312:   PetscCall(VecRestoreArray(xx, &x));
313:   PetscCall(VecRestoreArrayRead(bb, &b));
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }