1 /*$Id: aijfact.c,v 1.124 1999/10/24 14:02:14 bsmith Exp bsmith $*/ 2 3 #include "src/mat/impls/aij/seq/aij.h" 4 #include "src/vec/vecimpl.h" 5 #include "src/inline/dot.h" 6 7 #undef __FUNC__ 8 #define __FUNC__ "MatOrdering_Flow_SeqAIJ" 9 int MatOrdering_Flow_SeqAIJ(Mat mat,MatOrderingType type,IS *irow,IS *icol) 10 { 11 PetscFunctionBegin; 12 13 SETERRQ(PETSC_ERR_SUP,0,"Code not written"); 14 #if !defined(PETSC_USE_DEBUG) 15 PetscFunctionReturn(0); 16 #endif 17 } 18 19 /* 20 Factorization code for AIJ format. 21 */ 22 #undef __FUNC__ 23 #define __FUNC__ "MatLUFactorSymbolic_SeqAIJ" 24 int MatLUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,double f,Mat *B) 25 { 26 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b; 27 IS isicol; 28 int *r,*ic, ierr, i, n = a->m, *ai = a->i, *aj = a->j; 29 int *ainew,*ajnew, jmax,*fill, *ajtmp, nz,shift = a->indexshift; 30 int *idnew, idx, row,m,fm, nnz, nzi, realloc = 0,nzbd,*im; 31 32 PetscFunctionBegin; 33 PetscValidHeaderSpecific(isrow,IS_COOKIE); 34 PetscValidHeaderSpecific(iscol,IS_COOKIE); 35 if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square"); 36 37 ierr = ISInvertPermutation(iscol,&isicol);CHKERRQ(ierr); 38 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 39 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 40 41 /* get new row pointers */ 42 ainew = (int *) PetscMalloc( (n+1)*sizeof(int) );CHKPTRQ(ainew); 43 ainew[0] = -shift; 44 /* don't know how many column pointers are needed so estimate */ 45 jmax = (int) (f*ai[n]+(!shift)); 46 ajnew = (int *) PetscMalloc( (jmax)*sizeof(int) );CHKPTRQ(ajnew); 47 /* fill is a linked list of nonzeros in active row */ 48 fill = (int *) PetscMalloc( (2*n+1)*sizeof(int));CHKPTRQ(fill); 49 im = fill + n + 1; 50 /* idnew is location of diagonal in factor */ 51 idnew = (int *) PetscMalloc( (n+1)*sizeof(int));CHKPTRQ(idnew); 52 idnew[0] = -shift; 53 54 for ( i=0; i<n; i++ ) { 55 /* first copy previous fill into linked list */ 56 nnz = nz = ai[r[i]+1] - ai[r[i]]; 57 if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,1,"Empty row in matrix"); 58 ajtmp = aj + ai[r[i]] + shift; 59 fill[n] = n; 60 while (nz--) { 61 fm = n; 62 idx = ic[*ajtmp++ + shift]; 63 do { 64 m = fm; 65 fm = fill[m]; 66 } while (fm < idx); 67 fill[m] = idx; 68 fill[idx] = fm; 69 } 70 row = fill[n]; 71 while ( row < i ) { 72 ajtmp = ajnew + idnew[row] + (!shift); 73 nzbd = 1 + idnew[row] - ainew[row]; 74 nz = im[row] - nzbd; 75 fm = row; 76 while (nz-- > 0) { 77 idx = *ajtmp++ + shift; 78 nzbd++; 79 if (idx == i) im[row] = nzbd; 80 do { 81 m = fm; 82 fm = fill[m]; 83 } while (fm < idx); 84 if (fm != idx) { 85 fill[m] = idx; 86 fill[idx] = fm; 87 fm = idx; 88 nnz++; 89 } 90 } 91 row = fill[row]; 92 } 93 /* copy new filled row into permanent storage */ 94 ainew[i+1] = ainew[i] + nnz; 95 if (ainew[i+1] > jmax) { 96 97 /* estimate how much additional space we will need */ 98 /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */ 99 /* just double the memory each time */ 100 int maxadd = jmax; 101 /* maxadd = (int) ((f*(ai[n]+(!shift))*(n-i+5))/n); */ 102 if (maxadd < nnz) maxadd = (n-i)*(nnz+1); 103 jmax += maxadd; 104 105 /* allocate a longer ajnew */ 106 ajtmp = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(ajtmp); 107 ierr = PetscMemcpy(ajtmp,ajnew,(ainew[i]+shift)*sizeof(int));CHKERRQ(ierr); 108 ierr = PetscFree(ajnew);CHKERRQ(ierr); 109 ajnew = ajtmp; 110 realloc++; /* count how many times we realloc */ 111 } 112 ajtmp = ajnew + ainew[i] + shift; 113 fm = fill[n]; 114 nzi = 0; 115 im[i] = nnz; 116 while (nnz--) { 117 if (fm < i) nzi++; 118 *ajtmp++ = fm - shift; 119 fm = fill[fm]; 120 } 121 idnew[i] = ainew[i] + nzi; 122 } 123 if (ai[n] != 0) { 124 double af = ((double)ainew[n])/((double)ai[n]); 125 PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n", 126 realloc,f,af); 127 PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:Run with -pc_lu_fill %g or use \n",af); 128 PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:PCLUSetFill(pc,%g);\n",af); 129 PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ:for best performance.\n"); 130 } else { 131 PLogInfo(A,"MatLUFactorSymbolic_SeqAIJ: Empty matrix\n"); 132 } 133 134 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 135 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 136 137 ierr = PetscFree(fill);CHKERRQ(ierr); 138 139 /* put together the new matrix */ 140 ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,B);CHKERRQ(ierr); 141 PLogObjectParent(*B,isicol); 142 b = (Mat_SeqAIJ *) (*B)->data; 143 ierr = PetscFree(b->imax);CHKERRQ(ierr); 144 b->singlemalloc = 0; 145 /* the next line frees the default space generated by the Create() */ 146 ierr = PetscFree(b->a);CHKERRQ(ierr); 147 ierr = PetscFree(b->ilen);CHKERRQ(ierr); 148 b->a = (Scalar *) PetscMalloc((ainew[n]+shift+1)*sizeof(Scalar));CHKPTRQ(b->a); 149 b->j = ajnew; 150 b->i = ainew; 151 b->diag = idnew; 152 b->ilen = 0; 153 b->imax = 0; 154 b->row = isrow; 155 b->col = iscol; 156 b->icol = isicol; 157 b->solve_work = (Scalar *) PetscMalloc( (n+1)*sizeof(Scalar));CHKPTRQ(b->solve_work); 158 /* In b structure: Free imax, ilen, old a, old j. 159 Allocate idnew, solve_work, new a, new j */ 160 PLogObjectMemory(*B,(ainew[n]+shift-n)*(sizeof(int)+sizeof(Scalar))); 161 b->maxnz = b->nz = ainew[n] + shift; 162 163 (*B)->factor = FACTOR_LU;; 164 (*B)->info.factor_mallocs = realloc; 165 (*B)->info.fill_ratio_given = f; 166 (*B)->ops->lufactornumeric = A->ops->lufactornumeric; /* Use Inode variant if A has inodes */ 167 168 if (ai[n] != 0) { 169 (*B)->info.fill_ratio_needed = ((double)ainew[n])/((double)ai[n]); 170 } else { 171 (*B)->info.fill_ratio_needed = 0.0; 172 } 173 PetscFunctionReturn(0); 174 } 175 /* ----------------------------------------------------------- */ 176 extern int Mat_AIJ_CheckInode(Mat); 177 178 #undef __FUNC__ 179 #define __FUNC__ "MatLUFactorNumeric_SeqAIJ" 180 int MatLUFactorNumeric_SeqAIJ(Mat A,Mat *B) 181 { 182 Mat C = *B; 183 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b = (Mat_SeqAIJ *)C->data; 184 IS isrow = b->row, isicol = b->icol; 185 int *r,*ic, ierr, i, j, n = a->m, *ai = b->i, *aj = b->j; 186 int *ajtmpold, *ajtmp, nz, row, *ics, shift = a->indexshift; 187 int *diag_offset = b->diag,diag,k; 188 int preserve_row_sums = (int) a->ilu_preserve_row_sums; 189 register int *pj; 190 Scalar *rtmp,*v, *pc, multiplier,sum,inner_sum,*rowsums = 0; 191 double ssum; 192 register Scalar *pv, *rtmps,*u_values; 193 194 PetscFunctionBegin; 195 196 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 197 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 198 rtmp = (Scalar *) PetscMalloc( (n+1)*sizeof(Scalar) );CHKPTRQ(rtmp); 199 ierr = PetscMemzero(rtmp,(n+1)*sizeof(Scalar));CHKERRQ(ierr); 200 rtmps = rtmp + shift; ics = ic + shift; 201 202 /* precalculate row sums */ 203 if (preserve_row_sums) { 204 rowsums = (Scalar *) PetscMalloc( n*sizeof(Scalar) );CHKPTRQ(rowsums); 205 for ( i=0; i<n; i++ ) { 206 nz = a->i[r[i]+1] - a->i[r[i]]; 207 v = a->a + a->i[r[i]] + shift; 208 sum = 0.0; 209 for ( j=0; j<nz; j++ ) sum += v[j]; 210 rowsums[i] = sum; 211 } 212 } 213 214 for ( i=0; i<n; i++ ) { 215 nz = ai[i+1] - ai[i]; 216 ajtmp = aj + ai[i] + shift; 217 for ( j=0; j<nz; j++ ) rtmps[ajtmp[j]] = 0.0; 218 219 /* load in initial (unfactored row) */ 220 nz = a->i[r[i]+1] - a->i[r[i]]; 221 ajtmpold = a->j + a->i[r[i]] + shift; 222 v = a->a + a->i[r[i]] + shift; 223 for ( j=0; j<nz; j++ ) rtmp[ics[ajtmpold[j]]] = v[j]; 224 225 row = *ajtmp++ + shift; 226 while (row < i ) { 227 pc = rtmp + row; 228 if (*pc != 0.0) { 229 pv = b->a + diag_offset[row] + shift; 230 pj = b->j + diag_offset[row] + (!shift); 231 multiplier = *pc / *pv++; 232 *pc = multiplier; 233 nz = ai[row+1] - diag_offset[row] - 1; 234 for (j=0; j<nz; j++) rtmps[pj[j]] -= multiplier * pv[j]; 235 PLogFlops(2*nz); 236 } 237 row = *ajtmp++ + shift; 238 } 239 /* finished row so stick it into b->a */ 240 pv = b->a + ai[i] + shift; 241 pj = b->j + ai[i] + shift; 242 nz = ai[i+1] - ai[i]; 243 for ( j=0; j<nz; j++ ) {pv[j] = rtmps[pj[j]];} 244 diag = diag_offset[i] - ai[i]; 245 /* 246 Possibly adjust diagonal entry on current row to force 247 LU matrix to have same row sum as initial matrix. 248 */ 249 if (pv[diag] == 0.0) { 250 SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,0,"Zero pivot row %d",i); 251 } 252 if (preserve_row_sums) { 253 pj = b->j + ai[i] + shift; 254 sum = rowsums[i]; 255 for ( j=0; j<diag; j++ ) { 256 u_values = b->a + diag_offset[pj[j]] + shift; 257 nz = ai[pj[j]+1] - diag_offset[pj[j]]; 258 inner_sum = 0.0; 259 for ( k=0; k<nz; k++ ) { 260 inner_sum += u_values[k]; 261 } 262 sum -= pv[j]*inner_sum; 263 264 } 265 nz = ai[i+1] - diag_offset[i] - 1; 266 u_values = b->a + diag_offset[i] + 1 + shift; 267 for ( k=0; k<nz; k++ ) { 268 sum -= u_values[k]; 269 } 270 ssum = PetscAbsScalar(sum/pv[diag]); 271 if (ssum < 1000. && ssum > .001) pv[diag] = sum; 272 } 273 /* check pivot entry for current row */ 274 } 275 276 /* invert diagonal entries for simplier triangular solves */ 277 for ( i=0; i<n; i++ ) { 278 b->a[diag_offset[i]+shift] = 1.0/b->a[diag_offset[i]+shift]; 279 } 280 281 if (preserve_row_sums) {ierr = PetscFree(rowsums);CHKERRQ(ierr);} 282 ierr = PetscFree(rtmp);CHKERRQ(ierr); 283 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 284 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 285 C->factor = FACTOR_LU; 286 ierr = Mat_AIJ_CheckInode(C);CHKERRQ(ierr); 287 C->assembled = PETSC_TRUE; 288 PLogFlops(b->n); 289 PetscFunctionReturn(0); 290 } 291 /* ----------------------------------------------------------- */ 292 #undef __FUNC__ 293 #define __FUNC__ "MatLUFactor_SeqAIJ" 294 int MatLUFactor_SeqAIJ(Mat A,IS row,IS col,double f) 295 { 296 Mat_SeqAIJ *mat = (Mat_SeqAIJ *) A->data; 297 int ierr; 298 Mat C; 299 PetscOps *Abops; 300 MatOps Aops; 301 302 PetscFunctionBegin; 303 ierr = MatLUFactorSymbolic(A,row,col,f,&C);CHKERRQ(ierr); 304 ierr = MatLUFactorNumeric(A,&C);CHKERRQ(ierr); 305 306 /* free all the data structures from mat */ 307 ierr = PetscFree(mat->a);CHKERRQ(ierr); 308 if (!mat->singlemalloc) { 309 ierr = PetscFree(mat->i);CHKERRQ(ierr); 310 ierr = PetscFree(mat->j);CHKERRQ(ierr); 311 } 312 if (mat->diag) {ierr = PetscFree(mat->diag);CHKERRQ(ierr);} 313 if (mat->ilen) {ierr = PetscFree(mat->ilen);CHKERRQ(ierr);} 314 if (mat->imax) {ierr = PetscFree(mat->imax);CHKERRQ(ierr);} 315 if (mat->solve_work) {ierr = PetscFree(mat->solve_work);CHKERRQ(ierr);} 316 if (mat->inode.size) {ierr = PetscFree(mat->inode.size);CHKERRQ(ierr);} 317 if (mat->icol) {ierr = ISDestroy(mat->icol);CHKERRQ(ierr);} 318 ierr = PetscFree(mat);CHKERRQ(ierr); 319 320 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 321 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 322 323 /* 324 This is horrible, horrible code. We need to keep the 325 A pointers for the bops and ops but copy everything 326 else from C. 327 */ 328 Abops = A->bops; 329 Aops = A->ops; 330 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 331 mat = (Mat_SeqAIJ *) A->data; 332 PLogObjectParent(A,mat->icol); 333 334 A->bops = Abops; 335 A->ops = Aops; 336 A->qlist = 0; 337 338 PetscHeaderDestroy(C); 339 PetscFunctionReturn(0); 340 } 341 /* ----------------------------------------------------------- */ 342 #undef __FUNC__ 343 #define __FUNC__ "MatSolve_SeqAIJ" 344 int MatSolve_SeqAIJ(Mat A,Vec bb, Vec xx) 345 { 346 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 347 IS iscol = a->col, isrow = a->row; 348 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 349 int nz,shift = a->indexshift,*rout,*cout; 350 Scalar *x,*b,*tmp, *tmps, *aa = a->a, sum, *v; 351 352 PetscFunctionBegin; 353 if (!n) PetscFunctionReturn(0); 354 355 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 356 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 357 tmp = a->solve_work; 358 359 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 360 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1); 361 362 /* forward solve the lower triangular */ 363 tmp[0] = b[*r++]; 364 tmps = tmp + shift; 365 for ( i=1; i<n; i++ ) { 366 v = aa + ai[i] + shift; 367 vi = aj + ai[i] + shift; 368 nz = a->diag[i] - ai[i]; 369 sum = b[*r++]; 370 while (nz--) sum -= *v++ * tmps[*vi++]; 371 tmp[i] = sum; 372 } 373 374 /* backward solve the upper triangular */ 375 for ( i=n-1; i>=0; i-- ){ 376 v = aa + a->diag[i] + (!shift); 377 vi = aj + a->diag[i] + (!shift); 378 nz = ai[i+1] - a->diag[i] - 1; 379 sum = tmp[i]; 380 while (nz--) sum -= *v++ * tmps[*vi++]; 381 x[*c--] = tmp[i] = sum*aa[a->diag[i]+shift]; 382 } 383 384 ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr); 385 ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr); 386 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 387 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 388 PLogFlops(2*a->nz - a->n); 389 PetscFunctionReturn(0); 390 } 391 392 /* ----------------------------------------------------------- */ 393 #undef __FUNC__ 394 #define __FUNC__ "MatSolve_SeqAIJ_NaturalOrdering" 395 int MatSolve_SeqAIJ_NaturalOrdering(Mat A,Vec bb, Vec xx) 396 { 397 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 398 int n = a->m, *ai = a->i, *aj = a->j, *adiag = a->diag,ierr; 399 Scalar *x,*b, *aa = a->a, sum; 400 #if !defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ) 401 int adiag_i,i,*vi,nz,ai_i; 402 Scalar *v; 403 #endif 404 405 PetscFunctionBegin; 406 if (!n) PetscFunctionReturn(0); 407 if (a->indexshift) { 408 ierr = MatSolve_SeqAIJ(A,bb,xx);CHKERRQ(ierr); 409 PetscFunctionReturn(0); 410 } 411 412 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 413 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 414 415 #if defined(PETSC_USE_FORTRAN_KERNEL_SOLVEAIJ) 416 fortransolveaij_(&n,x,ai,aj,adiag,aa,b); 417 #else 418 /* forward solve the lower triangular */ 419 x[0] = b[0]; 420 for ( i=1; i<n; i++ ) { 421 ai_i = ai[i]; 422 v = aa + ai_i; 423 vi = aj + ai_i; 424 nz = adiag[i] - ai_i; 425 sum = b[i]; 426 while (nz--) sum -= *v++ * x[*vi++]; 427 x[i] = sum; 428 } 429 430 /* backward solve the upper triangular */ 431 for ( i=n-1; i>=0; i-- ){ 432 adiag_i = adiag[i]; 433 v = aa + adiag_i + 1; 434 vi = aj + adiag_i + 1; 435 nz = ai[i+1] - adiag_i - 1; 436 sum = x[i]; 437 while (nz--) sum -= *v++ * x[*vi++]; 438 x[i] = sum*aa[adiag_i]; 439 } 440 #endif 441 PLogFlops(2*a->nz - a->n); 442 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 443 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 444 PetscFunctionReturn(0); 445 } 446 447 #undef __FUNC__ 448 #define __FUNC__ "MatSolveAdd_SeqAIJ" 449 int MatSolveAdd_SeqAIJ(Mat A,Vec bb, Vec yy, Vec xx) 450 { 451 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 452 IS iscol = a->col, isrow = a->row; 453 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 454 int nz, shift = a->indexshift,*rout,*cout; 455 Scalar *x,*b,*tmp, *aa = a->a, sum, *v; 456 457 PetscFunctionBegin; 458 if (yy != xx) {ierr = VecCopy(yy,xx);CHKERRQ(ierr);} 459 460 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 461 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 462 tmp = a->solve_work; 463 464 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 465 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1); 466 467 /* forward solve the lower triangular */ 468 tmp[0] = b[*r++]; 469 for ( i=1; i<n; i++ ) { 470 v = aa + ai[i] + shift; 471 vi = aj + ai[i] + shift; 472 nz = a->diag[i] - ai[i]; 473 sum = b[*r++]; 474 while (nz--) sum -= *v++ * tmp[*vi++ + shift]; 475 tmp[i] = sum; 476 } 477 478 /* backward solve the upper triangular */ 479 for ( i=n-1; i>=0; i-- ){ 480 v = aa + a->diag[i] + (!shift); 481 vi = aj + a->diag[i] + (!shift); 482 nz = ai[i+1] - a->diag[i] - 1; 483 sum = tmp[i]; 484 while (nz--) sum -= *v++ * tmp[*vi++ + shift]; 485 tmp[i] = sum*aa[a->diag[i]+shift]; 486 x[*c--] += tmp[i]; 487 } 488 489 ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr); 490 ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr); 491 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 492 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 493 PLogFlops(2*a->nz); 494 495 PetscFunctionReturn(0); 496 } 497 /* -------------------------------------------------------------------*/ 498 #undef __FUNC__ 499 #define __FUNC__ "MatSolveTrans_SeqAIJ" 500 int MatSolveTrans_SeqAIJ(Mat A,Vec bb, Vec xx) 501 { 502 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 503 IS iscol = a->col, isrow = a->row; 504 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 505 int nz,shift = a->indexshift,*rout,*cout, *diag = a->diag; 506 Scalar *x,*b,*tmp, *aa = a->a, *v, s1; 507 508 PetscFunctionBegin; 509 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 510 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 511 tmp = a->solve_work; 512 513 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 514 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout; 515 516 /* copy the b into temp work space according to permutation */ 517 for ( i=0; i<n; i++ ) tmp[i] = b[c[i]]; 518 519 /* forward solve the U^T */ 520 for ( i=0; i<n; i++ ) { 521 v = aa + diag[i] + shift; 522 vi = aj + diag[i] + (!shift); 523 nz = ai[i+1] - diag[i] - 1; 524 s1 = tmp[i] *= *(v++); /* multiply by inverse of diagonal entry */ 525 while (nz--) { 526 tmp[*vi++ + shift] -= (*v++)*s1; 527 } 528 } 529 530 /* backward solve the L^T */ 531 for ( i=n-1; i>=0; i-- ){ 532 v = aa + diag[i] - 1 + shift; 533 vi = aj + diag[i] - 1 + shift; 534 nz = diag[i] - ai[i]; 535 s1 = tmp[i]; 536 while (nz--) { 537 tmp[*vi-- + shift] -= (*v--)*s1; 538 } 539 } 540 541 /* copy tmp into x according to permutation */ 542 for ( i=0; i<n; i++ ) x[r[i]] = tmp[i]; 543 544 ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr); 545 ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr); 546 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 547 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 548 549 PLogFlops(2*a->nz-a->n); 550 PetscFunctionReturn(0); 551 } 552 553 #undef __FUNC__ 554 #define __FUNC__ "MatSolveTransAdd_SeqAIJ" 555 int MatSolveTransAdd_SeqAIJ(Mat A,Vec bb, Vec zz,Vec xx) 556 { 557 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 558 IS iscol = a->col, isrow = a->row; 559 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 560 int nz,shift = a->indexshift, *rout, *cout, *diag = a->diag; 561 Scalar *x,*b,*tmp, *aa = a->a, *v; 562 563 PetscFunctionBegin; 564 if (zz != xx) VecCopy(zz,xx); 565 566 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 567 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 568 tmp = a->solve_work; 569 570 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 571 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout; 572 573 /* copy the b into temp work space according to permutation */ 574 for ( i=0; i<n; i++ ) tmp[i] = b[c[i]]; 575 576 /* forward solve the U^T */ 577 for ( i=0; i<n; i++ ) { 578 v = aa + diag[i] + shift; 579 vi = aj + diag[i] + (!shift); 580 nz = ai[i+1] - diag[i] - 1; 581 tmp[i] *= *v++; 582 while (nz--) { 583 tmp[*vi++ + shift] -= (*v++)*tmp[i]; 584 } 585 } 586 587 /* backward solve the L^T */ 588 for ( i=n-1; i>=0; i-- ){ 589 v = aa + diag[i] - 1 + shift; 590 vi = aj + diag[i] - 1 + shift; 591 nz = diag[i] - ai[i]; 592 while (nz--) { 593 tmp[*vi-- + shift] -= (*v--)*tmp[i]; 594 } 595 } 596 597 /* copy tmp into x according to permutation */ 598 for ( i=0; i<n; i++ ) x[r[i]] += tmp[i]; 599 600 ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr); 601 ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr); 602 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 603 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 604 605 PLogFlops(2*a->nz); 606 PetscFunctionReturn(0); 607 } 608 /* ----------------------------------------------------------------*/ 609 extern int MatMissingDiag_SeqAIJ(Mat); 610 611 #undef __FUNC__ 612 #define __FUNC__ "MatILUFactorSymbolic_SeqAIJ" 613 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,MatILUInfo *info,Mat *fact) 614 { 615 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b; 616 IS isicol; 617 int *r,*ic, ierr, prow, n = a->m, *ai = a->i, *aj = a->j; 618 int *ainew,*ajnew, jmax,*fill, *xi, nz, *im,*ajfill,*flev; 619 int *dloc, idx, row,m,fm, nzf, nzi,len, realloc = 0, dcount = 0; 620 int incrlev,nnz,i,shift = a->indexshift,levels,diagonal_fill; 621 PetscTruth col_identity, row_identity; 622 double f; 623 624 PetscFunctionBegin; 625 if (info) { 626 f = info->fill; 627 levels = (int) info->levels; 628 diagonal_fill = (int) info->diagonal_fill; 629 } else { 630 f = 1.0; 631 levels = 0; 632 diagonal_fill = 0; 633 } 634 ierr = ISInvertPermutation(iscol,&isicol);CHKERRQ(ierr); 635 636 /* special case that simply copies fill pattern */ 637 ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr); 638 ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr); 639 if (levels == 0 && row_identity && col_identity) { 640 ierr = MatDuplicate_SeqAIJ(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr); 641 (*fact)->factor = FACTOR_LU; 642 b = (Mat_SeqAIJ *) (*fact)->data; 643 if (!b->diag) { 644 ierr = MatMarkDiag_SeqAIJ(*fact);CHKERRQ(ierr); 645 } 646 ierr = MatMissingDiag_SeqAIJ(*fact);CHKERRQ(ierr); 647 b->row = isrow; 648 b->col = iscol; 649 b->icol = isicol; 650 b->solve_work = (Scalar *) PetscMalloc((b->m+1)*sizeof(Scalar));CHKPTRQ(b->solve_work); 651 (*fact)->ops->solve = MatSolve_SeqAIJ_NaturalOrdering; 652 PetscFunctionReturn(0); 653 } 654 655 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 656 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 657 658 /* get new row pointers */ 659 ainew = (int *) PetscMalloc( (n+1)*sizeof(int) );CHKPTRQ(ainew); 660 ainew[0] = -shift; 661 /* don't know how many column pointers are needed so estimate */ 662 jmax = (int) (f*(ai[n]+!shift)); 663 ajnew = (int *) PetscMalloc( (jmax)*sizeof(int) );CHKPTRQ(ajnew); 664 /* ajfill is level of fill for each fill entry */ 665 ajfill = (int *) PetscMalloc( (jmax)*sizeof(int) );CHKPTRQ(ajfill); 666 /* fill is a linked list of nonzeros in active row */ 667 fill = (int *) PetscMalloc( (n+1)*sizeof(int));CHKPTRQ(fill); 668 /* im is level for each filled value */ 669 im = (int *) PetscMalloc( (n+1)*sizeof(int));CHKPTRQ(im); 670 /* dloc is location of diagonal in factor */ 671 dloc = (int *) PetscMalloc( (n+1)*sizeof(int));CHKPTRQ(dloc); 672 dloc[0] = 0; 673 for ( prow=0; prow<n; prow++ ) { 674 675 /* copy current row into linked list */ 676 nzf = nz = ai[r[prow]+1] - ai[r[prow]]; 677 if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,1,"Empty row in matrix"); 678 xi = aj + ai[r[prow]] + shift; 679 fill[n] = n; 680 fill[prow] = -1; /* marker to indicate if diagonal exists */ 681 while (nz--) { 682 fm = n; 683 idx = ic[*xi++ + shift]; 684 do { 685 m = fm; 686 fm = fill[m]; 687 } while (fm < idx); 688 fill[m] = idx; 689 fill[idx] = fm; 690 im[idx] = 0; 691 } 692 693 /* make sure diagonal entry is included */ 694 if (diagonal_fill && fill[prow] == -1) { 695 fm = n; 696 while (fill[fm] < prow) fm = fill[fm]; 697 fill[prow] = fill[fm]; /* insert diagonal into linked list */ 698 fill[fm] = prow; 699 im[prow] = 0; 700 nzf++; 701 dcount++; 702 } 703 704 nzi = 0; 705 row = fill[n]; 706 while ( row < prow ) { 707 incrlev = im[row] + 1; 708 nz = dloc[row]; 709 xi = ajnew + ainew[row] + shift + nz + 1; 710 flev = ajfill + ainew[row] + shift + nz + 1; 711 nnz = ainew[row+1] - ainew[row] - nz - 1; 712 fm = row; 713 while (nnz-- > 0) { 714 idx = *xi++ + shift; 715 if (*flev + incrlev > levels) { 716 flev++; 717 continue; 718 } 719 do { 720 m = fm; 721 fm = fill[m]; 722 } while (fm < idx); 723 if (fm != idx) { 724 im[idx] = *flev + incrlev; 725 fill[m] = idx; 726 fill[idx] = fm; 727 fm = idx; 728 nzf++; 729 } else { 730 if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev; 731 } 732 flev++; 733 } 734 row = fill[row]; 735 nzi++; 736 } 737 /* copy new filled row into permanent storage */ 738 ainew[prow+1] = ainew[prow] + nzf; 739 if (ainew[prow+1] > jmax-shift) { 740 741 /* estimate how much additional space we will need */ 742 /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */ 743 /* just double the memory each time */ 744 /* maxadd = (int) ((f*(ai[n]+!shift)*(n-prow+5))/n); */ 745 int maxadd = jmax; 746 if (maxadd < nzf) maxadd = (n-prow)*(nzf+1); 747 jmax += maxadd; 748 749 /* allocate a longer ajnew and ajfill */ 750 xi = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(xi); 751 ierr = PetscMemcpy(xi,ajnew,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr); 752 ierr = PetscFree(ajnew);CHKERRQ(ierr); 753 ajnew = xi; 754 xi = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(xi); 755 ierr = PetscMemcpy(xi,ajfill,(ainew[prow]+shift)*sizeof(int));CHKERRQ(ierr); 756 ierr = PetscFree(ajfill);CHKERRQ(ierr); 757 ajfill = xi; 758 realloc++; /* count how many times we realloc */ 759 } 760 xi = ajnew + ainew[prow] + shift; 761 flev = ajfill + ainew[prow] + shift; 762 dloc[prow] = nzi; 763 fm = fill[n]; 764 while (nzf--) { 765 *xi++ = fm - shift; 766 *flev++ = im[fm]; 767 fm = fill[fm]; 768 } 769 /* make sure row has diagonal entry */ 770 if (ajnew[ainew[prow]+shift+dloc[prow]]+shift != prow) { 771 SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,1,"Row %d has missing diagonal in factored matrix\n\ 772 try running with -pc_ilu_nonzeros_along_diagonal or -pc_ilu_diagonal_fill",prow); 773 } 774 } 775 ierr = PetscFree(ajfill); CHKERRQ(ierr); 776 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 777 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 778 ierr = PetscFree(fill);CHKERRQ(ierr); 779 ierr = PetscFree(im);CHKERRQ(ierr); 780 781 { 782 double af = ((double)ainew[n])/((double)ai[n]); 783 PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n", 784 realloc,f,af); 785 PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Run with -pc_ilu_fill %g or use \n",af); 786 PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:PCILUSetFill(pc,%g);\n",af); 787 PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:for best performance.\n"); 788 if (diagonal_fill) { 789 PLogInfo(A,"MatILUFactorSymbolic_SeqAIJ:Detected and replace %d missing diagonals",dcount); 790 } 791 } 792 793 /* put together the new matrix */ 794 ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact);CHKERRQ(ierr); 795 PLogObjectParent(*fact,isicol); 796 b = (Mat_SeqAIJ *) (*fact)->data; 797 ierr = PetscFree(b->imax);CHKERRQ(ierr); 798 b->singlemalloc = 0; 799 len = (ainew[n] + shift)*sizeof(Scalar); 800 /* the next line frees the default space generated by the Create() */ 801 ierr = PetscFree(b->a);CHKERRQ(ierr); 802 ierr = PetscFree(b->ilen);CHKERRQ(ierr); 803 b->a = (Scalar *) PetscMalloc( len+1 );CHKPTRQ(b->a); 804 b->j = ajnew; 805 b->i = ainew; 806 for ( i=0; i<n; i++ ) dloc[i] += ainew[i]; 807 b->diag = dloc; 808 b->ilen = 0; 809 b->imax = 0; 810 b->row = isrow; 811 b->col = iscol; 812 b->icol = isicol; 813 b->solve_work = (Scalar *) PetscMalloc( (n+1)*sizeof(Scalar));CHKPTRQ(b->solve_work); 814 /* In b structure: Free imax, ilen, old a, old j. 815 Allocate dloc, solve_work, new a, new j */ 816 PLogObjectMemory(*fact,(ainew[n]+shift-n) * (sizeof(int)+sizeof(Scalar))); 817 b->maxnz = b->nz = ainew[n] + shift; 818 (*fact)->factor = FACTOR_LU; 819 820 (*fact)->info.factor_mallocs = realloc; 821 (*fact)->info.fill_ratio_given = f; 822 (*fact)->info.fill_ratio_needed = ((double)ainew[n])/((double)ai[prow]); 823 (*fact)->factor = FACTOR_LU;; 824 825 PetscFunctionReturn(0); 826 } 827 828 829 830 831