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