Actual source code: relax.h

petsc-3.11.3 2019-06-26
Report Typos and Errors

  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_Hermitian_ushort(Mat A,Vec xx,Vec zz)
 43: #else
 44: PetscErrorCode MatMult_SeqSBAIJ_1_Hermitian(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;

 65:   VecSet(zz,0.0);
 66:   VecGetArrayRead(xx,&x);
 67:   VecGetArray(zz,&z);

 69:   v = a->a;
 70:   for (i=0; i<mbs; i++) {
 71:     nz = ai[i+1] - ai[i];    /* length of i_th row of A */
 72:     if (!nz) continue; /* Move to the next row if the current row is empty */
 73:     nonzerorow++;
 74:     x1   = x[i];
 75:     sum  = 0.0;
 76:     jmin = 0;
 77:     if (ib[0] == i) {
 78:       sum = v[0]*x1;           /* diagonal term */
 79:       jmin++;
 80:     }
 81:     for (j=jmin; j<nz; j++) {
 82:       ibt     = ib[j];
 83:       vj      = v[j];
 84:       sum    += vj * x[ibt]; /* (strict upper triangular part of A)*x  */
 85:       z[ibt] += PetscConj(v[j]) * x1;    /* (strict lower triangular part of A)*x  */
 86:     }
 87:     z[i] += sum;
 88:     v    +=    nz;
 89:     ib   += nz;
 90:   }

 92:   VecRestoreArrayRead(xx,&x);
 93:   VecRestoreArray(zz,&z);
 94:   PetscLogFlops(2.0*(2.0*a->nz - nonzerorow) - nonzerorow);
 95:   return(0);
 96: }

 98: #if defined(USESHORT)
 99: PetscErrorCode MatMult_SeqSBAIJ_1_ushort(Mat A,Vec xx,Vec zz)
100: #else
101: PetscErrorCode MatMult_SeqSBAIJ_1(Mat A,Vec xx,Vec zz)
102: #endif
103: {
104:   Mat_SeqSBAIJ      *a = (Mat_SeqSBAIJ*)A->data;
105:   const PetscScalar *x;
106:   PetscScalar       *z,x1,sum;
107:   const MatScalar   *v;
108:   MatScalar         vj;
109:   PetscErrorCode    ierr;
110:   PetscInt          mbs=a->mbs,i,j,nz;
111:   const PetscInt    *ai=a->i;
112: #if defined(USESHORT)
113:   const unsigned short *ib=a->jshort;
114:   unsigned short       ibt;
115: #else
116:   const PetscInt *ib=a->j;
117:   PetscInt       ibt;
118: #endif
119:   PetscInt nonzerorow=0,jmin;

122:   VecSet(zz,0.0);
123:   VecGetArrayRead(xx,&x);
124:   VecGetArray(zz,&z);

126:   v = a->a;
127:   for (i=0; i<mbs; i++) {
128:     nz = ai[i+1] - ai[i];          /* length of i_th row of A */
129:     if (!nz) continue; /* Move to the next row if the current row is empty */
130:     nonzerorow++;
131:     sum  = 0.0;
132:     jmin = 0;
133:     x1   = x[i];
134:     if (ib[0] == i) {
135:       sum = v[0]*x1;                 /* diagonal term */
136:       jmin++;
137:     }
138:     PetscPrefetchBlock(ib+nz,nz,0,PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
139:     PetscPrefetchBlock(v+nz,nz,0,PETSC_PREFETCH_HINT_NTA);  /* Entries for the next row */
140:     for (j=jmin; j<nz; j++) {
141:       ibt     = ib[j];
142:       vj      = v[j];
143:       z[ibt] += vj * x1;       /* (strict lower triangular part of A)*x  */
144:       sum    += vj * x[ibt];   /* (strict upper triangular part of A)*x  */
145:     }
146:     z[i] += sum;
147:     v    += nz;
148:     ib   += nz;
149:   }

151:   VecRestoreArrayRead(xx,&x);
152:   VecRestoreArray(zz,&z);
153:   PetscLogFlops(2.0*(2.0*a->nz - nonzerorow) - nonzerorow);
154:   return(0);
155: }

157: #if defined(USESHORT)
158: PetscErrorCode MatSOR_SeqSBAIJ_ushort(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
159: #else
160: PetscErrorCode MatSOR_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
161: #endif
162: {
163:   Mat_SeqSBAIJ      *a = (Mat_SeqSBAIJ*)A->data;
164:   const MatScalar   *aa=a->a,*v,*v1,*aidiag;
165:   PetscScalar       *x,*t,sum;
166:   const PetscScalar *b;
167:   MatScalar         tmp;
168:   PetscErrorCode    ierr;
169:   PetscInt          m  =a->mbs,bs=A->rmap->bs,j;
170:   const PetscInt    *ai=a->i;
171: #if defined(USESHORT)
172:   const unsigned short *aj=a->jshort,*vj,*vj1;
173: #else
174:   const PetscInt *aj=a->j,*vj,*vj1;
175: #endif
176:   PetscInt nz,nz1,i;

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

182:   its = its*lits;
183:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);

185:   if (bs > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented");

187:   VecGetArray(xx,&x);
188:   VecGetArrayRead(bb,&b);

190:   if (!a->idiagvalid) {
191:     if (!a->idiag) {
192:       PetscMalloc1(m,&a->idiag);
193:     }
194:     for (i=0; i<a->mbs; i++) a->idiag[i] = 1.0/a->a[a->i[i]];
195:     a->idiagvalid = PETSC_TRUE;
196:   }

198:   if (!a->sor_work) {
199:     PetscMalloc1(m,&a->sor_work);
200:   }
201:   t = a->sor_work;

203:   aidiag = a->idiag;

205:   if (flag == SOR_APPLY_UPPER) {
206:     /* apply (U + D/omega) to the vector */
207:     PetscScalar d;
208:     for (i=0; i<m; i++) {
209:       d   = fshift + aa[ai[i]];
210:       nz  = ai[i+1] - ai[i] - 1;
211:       vj  = aj + ai[i] + 1;
212:       v   = aa + ai[i] + 1;
213:       sum = b[i]*d/omega;
214: #ifdef USESHORT
215:       PetscSparseDensePlusDot_no_function(sum,b,v,vj,nz);
216: #else
217:       PetscSparseDensePlusDot(sum,b,v,vj,nz);
218: #endif
219:       x[i] = sum;
220:     }
221:     PetscLogFlops(a->nz);
222:   }

224:   if (flag & SOR_ZERO_INITIAL_GUESS) {
225:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
226:       PetscMemcpy(t,b,m*sizeof(PetscScalar));

228:       v  = aa + 1;
229:       vj = aj + 1;
230:       for (i=0; i<m; i++) {
231:         nz  = ai[i+1] - ai[i] - 1;
232:         tmp = -(x[i] = omega*t[i]*aidiag[i]);
233:         for (j=0; j<nz; j++) t[vj[j]] += tmp*v[j];
234:         v  += nz + 1;
235:         vj += nz + 1;
236:       }
237:       PetscLogFlops(2*a->nz);
238:     }

240:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
241:       int nz2;
242:       if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)) {
243: #if defined(PETSC_USE_BACKWARD_LOOP)
244:         v  = aa + ai[m] - 1;
245:         vj = aj + ai[m] - 1;
246:         for (i=m-1; i>=0; i--) {
247:           sum = b[i];
248:           nz  = ai[i+1] - ai[i] - 1;
249:           {PetscInt __i;for (__i=0; __i<nz; __i++) sum -= v[-__i] * x[vj[-__i]];}
250: #else
251:         v  = aa + ai[m-1] + 1;
252:         vj = aj + ai[m-1] + 1;
253:         nz = 0;
254:         for (i=m-1; i>=0; i--) {
255:           sum = b[i];
256:           nz2 = ai[i] - ai[PetscMax(i-1,0)] - 1; /* avoid referencing ai[-1], nonsense nz2 is okay on last iteration */
257:           PETSC_Prefetch(v-nz2-1,0,PETSC_PREFETCH_HINT_NTA);
258:           PETSC_Prefetch(vj-nz2-1,0,PETSC_PREFETCH_HINT_NTA);
259:           PetscSparseDenseMinusDot(sum,x,v,vj,nz);
260:           nz = nz2;
261: #endif
262:           x[i] = omega*sum*aidiag[i];
263:           v   -= nz + 1;
264:           vj  -= nz + 1;
265:         }
266:         PetscLogFlops(2*a->nz);
267:       } else {
268:         v  = aa + ai[m-1] + 1;
269:         vj = aj + ai[m-1] + 1;
270:         nz = 0;
271:         for (i=m-1; i>=0; i--) {
272:           sum = t[i];
273:           nz2 = ai[i] - ai[PetscMax(i-1,0)] - 1; /* avoid referencing ai[-1], nonsense nz2 is okay on last iteration */
274:           PETSC_Prefetch(v-nz2-1,0,PETSC_PREFETCH_HINT_NTA);
275:           PETSC_Prefetch(vj-nz2-1,0,PETSC_PREFETCH_HINT_NTA);
276:           PetscSparseDenseMinusDot(sum,x,v,vj,nz);
277:           x[i] = (1-omega)*x[i] + omega*sum*aidiag[i];
278:           nz   = nz2;
279:           v   -= nz + 1;
280:           vj  -= nz + 1;
281:         }
282:         PetscLogFlops(2*a->nz);
283:       }
284:     }
285:     its--;
286:   }

288:   while (its--) {
289:     /*
290:        forward sweep:
291:        for i=0,...,m-1:
292:          sum[i] = (b[i] - U(i,:)x)/d[i];
293:          x[i]   = (1-omega)x[i] + omega*sum[i];
294:          b      = b - x[i]*U^T(i,:);

296:     */
297:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
298:       PetscMemcpy(t,b,m*sizeof(PetscScalar));

300:       for (i=0; i<m; i++) {
301:         v    = aa + ai[i] + 1; v1=v;
302:         vj   = aj + ai[i] + 1; vj1=vj;
303:         nz   = ai[i+1] - ai[i] - 1; nz1=nz;
304:         sum  = t[i];
305:         while (nz1--) sum -= (*v1++)*x[*vj1++];
306:         x[i] = (1-omega)*x[i] + omega*sum*aidiag[i];
307:         while (nz--) t[*vj++] -= x[i]*(*v++);
308:       }
309:       PetscLogFlops(4.0*a->nz);
310:     }

312:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
313:       /*
314:        backward sweep:
315:        b = b - x[i]*U^T(i,:), i=0,...,n-2
316:        for i=m-1,...,0:
317:          sum[i] = (b[i] - U(i,:)x)/d[i];
318:          x[i]   = (1-omega)x[i] + omega*sum[i];
319:       */
320:       /* if there was a forward sweep done above then I thing the next two for loops are not needed */
321:       PetscMemcpy(t,b,m*sizeof(PetscScalar));

323:       for (i=0; i<m-1; i++) {  /* update rhs */
324:         v    = aa + ai[i] + 1;
325:         vj   = aj + ai[i] + 1;
326:         nz   = ai[i+1] - ai[i] - 1;
327:         while (nz--) t[*vj++] -= x[i]*(*v++);
328:       }
329:       PetscLogFlops(2.0*(a->nz - m));
330:       for (i=m-1; i>=0; i--) {
331:         v    = aa + ai[i] + 1;
332:         vj   = aj + ai[i] + 1;
333:         nz   = ai[i+1] - ai[i] - 1;
334:         sum  = t[i];
335:         while (nz--) sum -= x[*vj++]*(*v++);
336:         x[i] =   (1-omega)*x[i] + omega*sum*aidiag[i];
337:       }
338:       PetscLogFlops(2.0*(a->nz + m));
339:     }
340:   }

342:   VecRestoreArray(xx,&x);
343:   VecRestoreArrayRead(bb,&b);
344:   return(0);
345: }