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) 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) 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