1 2 /* 3 Factorization code for BAIJ format. 4 */ 5 #include <../src/mat/impls/baij/seq/baij.h> 6 #include <../src/mat/blockinvert.h> 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "MatLUFactorNumeric_SeqBAIJ_2" 10 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2(Mat B,Mat A,const MatFactorInfo *info) 11 { 12 Mat C=B; 13 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ *)C->data; 14 IS isrow = b->row,isicol = b->icol; 15 PetscErrorCode ierr; 16 const PetscInt *r,*ic; 17 PetscInt i,j,k,nz,nzL,row,*pj; 18 const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2; 19 const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag; 20 MatScalar *rtmp,*pc,*mwork,*pv; 21 MatScalar *aa=a->a,*v; 22 PetscInt flg; 23 PetscReal shift = info->shiftamount; 24 25 PetscFunctionBegin; 26 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 27 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 28 29 /* generate work space needed by the factorization */ 30 ierr = PetscMalloc2(bs2*n,MatScalar,&rtmp,bs2,MatScalar,&mwork);CHKERRQ(ierr); 31 ierr = PetscMemzero(rtmp,bs2*n*sizeof(MatScalar));CHKERRQ(ierr); 32 33 for (i=0; i<n; i++) { 34 /* zero rtmp */ 35 /* L part */ 36 nz = bi[i+1] - bi[i]; 37 bjtmp = bj + bi[i]; 38 for (j=0; j<nz; j++) { 39 ierr = PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 40 } 41 42 /* U part */ 43 nz = bdiag[i] - bdiag[i+1]; 44 bjtmp = bj + bdiag[i+1]+1; 45 for (j=0; j<nz; j++) { 46 ierr = PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 47 } 48 49 /* load in initial (unfactored row) */ 50 nz = ai[r[i]+1] - ai[r[i]]; 51 ajtmp = aj + ai[r[i]]; 52 v = aa + bs2*ai[r[i]]; 53 for (j=0; j<nz; j++) { 54 ierr = PetscMemcpy(rtmp+bs2*ic[ajtmp[j]],v+bs2*j,bs2*sizeof(MatScalar));CHKERRQ(ierr); 55 } 56 57 /* elimination */ 58 bjtmp = bj + bi[i]; 59 nzL = bi[i+1] - bi[i]; 60 for (k=0;k < nzL;k++) { 61 row = bjtmp[k]; 62 pc = rtmp + bs2*row; 63 for (flg=0,j=0; j<bs2; j++) { 64 if (pc[j] != (PetscScalar)0.0) { 65 flg = 1; 66 break; 67 } 68 } 69 if (flg) { 70 pv = b->a + bs2*bdiag[row]; 71 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ 72 ierr = PetscKernel_A_gets_A_times_B_2(pc,pv,mwork);CHKERRQ(ierr); 73 74 pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */ 75 pv = b->a + bs2*(bdiag[row+1]+1); 76 nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */ 77 for (j=0; j<nz; j++) { 78 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ 79 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ 80 v = rtmp + 4*pj[j]; 81 ierr = PetscKernel_A_gets_A_minus_B_times_C_2(v,pc,pv);CHKERRQ(ierr); 82 pv += 4; 83 } 84 ierr = PetscLogFlops(16*nz+12);CHKERRQ(ierr); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 85 } 86 } 87 88 /* finished row so stick it into b->a */ 89 /* L part */ 90 pv = b->a + bs2*bi[i] ; 91 pj = b->j + bi[i] ; 92 nz = bi[i+1] - bi[i]; 93 for (j=0; j<nz; j++) { 94 ierr = PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 95 } 96 97 /* Mark diagonal and invert diagonal for simplier triangular solves */ 98 pv = b->a + bs2*bdiag[i]; 99 pj = b->j + bdiag[i]; 100 ierr = PetscMemcpy(pv,rtmp+bs2*pj[0],bs2*sizeof(MatScalar));CHKERRQ(ierr); 101 /* ierr = PetscKernel_A_gets_inverse_A(bs,pv,v_pivots,v_work);CHKERRQ(ierr); */ 102 ierr = PetscKernel_A_gets_inverse_A_2(pv,shift);CHKERRQ(ierr); 103 104 /* U part */ 105 pv = b->a + bs2*(bdiag[i+1]+1); 106 pj = b->j + bdiag[i+1]+1; 107 nz = bdiag[i] - bdiag[i+1] - 1; 108 for (j=0; j<nz; j++) { 109 ierr = PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 110 } 111 } 112 113 ierr = PetscFree2(rtmp,mwork);CHKERRQ(ierr); 114 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 115 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 116 C->ops->solve = MatSolve_SeqBAIJ_2; 117 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2; 118 119 C->assembled = PETSC_TRUE; 120 ierr = PetscLogFlops(1.333333333333*2*2*2*n);CHKERRQ(ierr); /* from inverting diagonal blocks */ 121 PetscFunctionReturn(0); 122 } 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering" 126 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info) 127 { 128 Mat C=B; 129 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ *)C->data; 130 PetscErrorCode ierr; 131 PetscInt i,j,k,nz,nzL,row,*pj; 132 const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2; 133 const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag; 134 MatScalar *rtmp,*pc,*mwork,*pv; 135 MatScalar *aa=a->a,*v; 136 PetscInt flg; 137 PetscReal shift = info->shiftamount; 138 139 PetscFunctionBegin; 140 /* generate work space needed by the factorization */ 141 ierr = PetscMalloc2(bs2*n,MatScalar,&rtmp,bs2,MatScalar,&mwork);CHKERRQ(ierr); 142 ierr = PetscMemzero(rtmp,bs2*n*sizeof(MatScalar));CHKERRQ(ierr); 143 144 for (i=0; i<n; i++) { 145 /* zero rtmp */ 146 /* L part */ 147 nz = bi[i+1] - bi[i]; 148 bjtmp = bj + bi[i]; 149 for (j=0; j<nz; j++) { 150 ierr = PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 151 } 152 153 /* U part */ 154 nz = bdiag[i] - bdiag[i+1]; 155 bjtmp = bj + bdiag[i+1]+1; 156 for (j=0; j<nz; j++) { 157 ierr = PetscMemzero(rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 158 } 159 160 /* load in initial (unfactored row) */ 161 nz = ai[i+1] - ai[i]; 162 ajtmp = aj + ai[i]; 163 v = aa + bs2*ai[i]; 164 for (j=0; j<nz; j++) { 165 ierr = PetscMemcpy(rtmp+bs2*ajtmp[j],v+bs2*j,bs2*sizeof(MatScalar));CHKERRQ(ierr); 166 } 167 168 /* elimination */ 169 bjtmp = bj + bi[i]; 170 nzL = bi[i+1] - bi[i]; 171 for (k=0;k < nzL;k++) { 172 row = bjtmp[k]; 173 pc = rtmp + bs2*row; 174 for (flg=0,j=0; j<bs2; j++) { 175 if (pc[j]!=(PetscScalar)0.0) { 176 flg = 1; 177 break; 178 } 179 } 180 if (flg) { 181 pv = b->a + bs2*bdiag[row]; 182 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ 183 ierr = PetscKernel_A_gets_A_times_B_2(pc,pv,mwork);CHKERRQ(ierr); 184 185 pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */ 186 pv = b->a + bs2*(bdiag[row+1]+1); 187 nz = bdiag[row]-bdiag[row+1] - 1; /* num of entries in U(row,:) excluding diag */ 188 for (j=0; j<nz; j++) { 189 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ 190 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ 191 v = rtmp + 4*pj[j]; 192 ierr = PetscKernel_A_gets_A_minus_B_times_C_2(v,pc,pv);CHKERRQ(ierr); 193 pv += 4; 194 } 195 ierr = PetscLogFlops(16*nz+12);CHKERRQ(ierr); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 196 } 197 } 198 199 /* finished row so stick it into b->a */ 200 /* L part */ 201 pv = b->a + bs2*bi[i] ; 202 pj = b->j + bi[i] ; 203 nz = bi[i+1] - bi[i]; 204 for (j=0; j<nz; j++) { 205 ierr = PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 206 } 207 208 /* Mark diagonal and invert diagonal for simplier triangular solves */ 209 pv = b->a + bs2*bdiag[i]; 210 pj = b->j + bdiag[i]; 211 ierr = PetscMemcpy(pv,rtmp+bs2*pj[0],bs2*sizeof(MatScalar));CHKERRQ(ierr); 212 /* ierr = PetscKernel_A_gets_inverse_A(bs,pv,v_pivots,v_work);CHKERRQ(ierr); */ 213 ierr = PetscKernel_A_gets_inverse_A_2(pv,shift);CHKERRQ(ierr); 214 215 /* U part */ 216 /* 217 pv = b->a + bs2*bi[2*n-i]; 218 pj = b->j + bi[2*n-i]; 219 nz = bi[2*n-i+1] - bi[2*n-i] - 1; 220 */ 221 pv = b->a + bs2*(bdiag[i+1]+1); 222 pj = b->j + bdiag[i+1]+1; 223 nz = bdiag[i] - bdiag[i+1] - 1; 224 for (j=0; j<nz; j++) { 225 ierr = PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 226 } 227 } 228 ierr = PetscFree2(rtmp,mwork);CHKERRQ(ierr); 229 230 C->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering; 231 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 232 C->assembled = PETSC_TRUE; 233 ierr = PetscLogFlops(1.333333333333*2*2*2*n);CHKERRQ(ierr); /* from inverting diagonal blocks */ 234 PetscFunctionReturn(0); 235 } 236 237 #undef __FUNCT__ 238 #define __FUNCT__ "MatLUFactorNumeric_SeqBAIJ_2_inplace" 239 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_inplace(Mat B,Mat A,const MatFactorInfo *info) 240 { 241 Mat C = B; 242 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ *)C->data; 243 IS isrow = b->row,isicol = b->icol; 244 PetscErrorCode ierr; 245 const PetscInt *r,*ic; 246 PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j; 247 PetscInt *ajtmpold,*ajtmp,nz,row; 248 PetscInt *diag_offset=b->diag,idx,*ai=a->i,*aj=a->j,*pj; 249 MatScalar *pv,*v,*rtmp,m1,m2,m3,m4,*pc,*w,*x,x1,x2,x3,x4; 250 MatScalar p1,p2,p3,p4; 251 MatScalar *ba = b->a,*aa = a->a; 252 PetscReal shift = info->shiftamount; 253 254 PetscFunctionBegin; 255 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 256 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 257 ierr = PetscMalloc(4*(n+1)*sizeof(MatScalar),&rtmp);CHKERRQ(ierr); 258 259 for (i=0; i<n; i++) { 260 nz = bi[i+1] - bi[i]; 261 ajtmp = bj + bi[i]; 262 for (j=0; j<nz; j++) { 263 x = rtmp+4*ajtmp[j]; x[0] = x[1] = x[2] = x[3] = 0.0; 264 } 265 /* load in initial (unfactored row) */ 266 idx = r[i]; 267 nz = ai[idx+1] - ai[idx]; 268 ajtmpold = aj + ai[idx]; 269 v = aa + 4*ai[idx]; 270 for (j=0; j<nz; j++) { 271 x = rtmp+4*ic[ajtmpold[j]]; 272 x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3]; 273 v += 4; 274 } 275 row = *ajtmp++; 276 while (row < i) { 277 pc = rtmp + 4*row; 278 p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3]; 279 if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) { 280 pv = ba + 4*diag_offset[row]; 281 pj = bj + diag_offset[row] + 1; 282 x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; 283 pc[0] = m1 = p1*x1 + p3*x2; 284 pc[1] = m2 = p2*x1 + p4*x2; 285 pc[2] = m3 = p1*x3 + p3*x4; 286 pc[3] = m4 = p2*x3 + p4*x4; 287 nz = bi[row+1] - diag_offset[row] - 1; 288 pv += 4; 289 for (j=0; j<nz; j++) { 290 x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; 291 x = rtmp + 4*pj[j]; 292 x[0] -= m1*x1 + m3*x2; 293 x[1] -= m2*x1 + m4*x2; 294 x[2] -= m1*x3 + m3*x4; 295 x[3] -= m2*x3 + m4*x4; 296 pv += 4; 297 } 298 ierr = PetscLogFlops(16.0*nz+12.0);CHKERRQ(ierr); 299 } 300 row = *ajtmp++; 301 } 302 /* finished row so stick it into b->a */ 303 pv = ba + 4*bi[i]; 304 pj = bj + bi[i]; 305 nz = bi[i+1] - bi[i]; 306 for (j=0; j<nz; j++) { 307 x = rtmp+4*pj[j]; 308 pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3]; 309 pv += 4; 310 } 311 /* invert diagonal block */ 312 w = ba + 4*diag_offset[i]; 313 ierr = PetscKernel_A_gets_inverse_A_2(w,shift);CHKERRQ(ierr); 314 } 315 316 ierr = PetscFree(rtmp);CHKERRQ(ierr); 317 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 318 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 319 C->ops->solve = MatSolve_SeqBAIJ_2_inplace; 320 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_inplace; 321 C->assembled = PETSC_TRUE; 322 ierr = PetscLogFlops(1.333333333333*8*b->mbs);CHKERRQ(ierr); /* from inverting diagonal blocks */ 323 PetscFunctionReturn(0); 324 } 325 /* 326 Version for when blocks are 2 by 2 Using natural ordering 327 */ 328 #undef __FUNCT__ 329 #define __FUNCT__ "MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace" 330 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info) 331 { 332 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ *)C->data; 333 PetscErrorCode ierr; 334 PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j; 335 PetscInt *ajtmpold,*ajtmp,nz,row; 336 PetscInt *diag_offset = b->diag,*ai=a->i,*aj=a->j,*pj; 337 MatScalar *pv,*v,*rtmp,*pc,*w,*x; 338 MatScalar p1,p2,p3,p4,m1,m2,m3,m4,x1,x2,x3,x4; 339 MatScalar *ba = b->a,*aa = a->a; 340 PetscReal shift = info->shiftamount; 341 342 PetscFunctionBegin; 343 ierr = PetscMalloc(4*(n+1)*sizeof(MatScalar),&rtmp);CHKERRQ(ierr); 344 for (i=0; i<n; i++) { 345 nz = bi[i+1] - bi[i]; 346 ajtmp = bj + bi[i]; 347 for (j=0; j<nz; j++) { 348 x = rtmp+4*ajtmp[j]; 349 x[0] = x[1] = x[2] = x[3] = 0.0; 350 } 351 /* load in initial (unfactored row) */ 352 nz = ai[i+1] - ai[i]; 353 ajtmpold = aj + ai[i]; 354 v = aa + 4*ai[i]; 355 for (j=0; j<nz; j++) { 356 x = rtmp+4*ajtmpold[j]; 357 x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3]; 358 v += 4; 359 } 360 row = *ajtmp++; 361 while (row < i) { 362 pc = rtmp + 4*row; 363 p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3]; 364 if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) { 365 pv = ba + 4*diag_offset[row]; 366 pj = bj + diag_offset[row] + 1; 367 x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; 368 pc[0] = m1 = p1*x1 + p3*x2; 369 pc[1] = m2 = p2*x1 + p4*x2; 370 pc[2] = m3 = p1*x3 + p3*x4; 371 pc[3] = m4 = p2*x3 + p4*x4; 372 nz = bi[row+1] - diag_offset[row] - 1; 373 pv += 4; 374 for (j=0; j<nz; j++) { 375 x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; 376 x = rtmp + 4*pj[j]; 377 x[0] -= m1*x1 + m3*x2; 378 x[1] -= m2*x1 + m4*x2; 379 x[2] -= m1*x3 + m3*x4; 380 x[3] -= m2*x3 + m4*x4; 381 pv += 4; 382 } 383 ierr = PetscLogFlops(16.0*nz+12.0);CHKERRQ(ierr); 384 } 385 row = *ajtmp++; 386 } 387 /* finished row so stick it into b->a */ 388 pv = ba + 4*bi[i]; 389 pj = bj + bi[i]; 390 nz = bi[i+1] - bi[i]; 391 for (j=0; j<nz; j++) { 392 x = rtmp+4*pj[j]; 393 pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3]; 394 /* 395 printf(" col %d:",pj[j]); 396 PetscInt j1; 397 for (j1=0; j1<4; j1++) printf(" %g,",*(pv+j1)); 398 printf("\n"); 399 */ 400 pv += 4; 401 } 402 /* invert diagonal block */ 403 w = ba + 4*diag_offset[i]; 404 /* 405 printf(" \n%d -th: diag: ",i); 406 for (j=0; j<4; j++) { 407 printf(" %g,",w[j]); 408 } 409 printf("\n----------------------------\n"); 410 */ 411 ierr = PetscKernel_A_gets_inverse_A_2(w,shift);CHKERRQ(ierr); 412 } 413 414 ierr = PetscFree(rtmp);CHKERRQ(ierr); 415 C->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering_inplace; 416 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace; 417 C->assembled = PETSC_TRUE; 418 ierr = PetscLogFlops(1.333333333333*8*b->mbs);CHKERRQ(ierr); /* from inverting diagonal blocks */ 419 PetscFunctionReturn(0); 420 } 421 422 /* ----------------------------------------------------------- */ 423 /* 424 Version for when blocks are 1 by 1. 425 */ 426 #undef __FUNCT__ 427 #define __FUNCT__ "MatLUFactorNumeric_SeqBAIJ_1" 428 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1(Mat B,Mat A,const MatFactorInfo *info) 429 { 430 Mat C=B; 431 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ *)C->data; 432 IS isrow = b->row,isicol = b->icol; 433 PetscErrorCode ierr; 434 const PetscInt *r,*ic,*ics; 435 const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bdiag=b->diag; 436 PetscInt i,j,k,nz,nzL,row,*pj; 437 const PetscInt *ajtmp,*bjtmp; 438 MatScalar *rtmp,*pc,multiplier,*pv; 439 const MatScalar *aa=a->a,*v; 440 PetscBool row_identity,col_identity; 441 FactorShiftCtx sctx; 442 const PetscInt *ddiag; 443 PetscReal rs; 444 MatScalar d; 445 446 PetscFunctionBegin; 447 /* MatPivotSetUp(): initialize shift context sctx */ 448 ierr = PetscMemzero(&sctx,sizeof(FactorShiftCtx));CHKERRQ(ierr); 449 450 if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */ 451 ddiag = a->diag; 452 sctx.shift_top = info->zeropivot; 453 for (i=0; i<n; i++) { 454 /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */ 455 d = (aa)[ddiag[i]]; 456 rs = -PetscAbsScalar(d) - PetscRealPart(d); 457 v = aa+ai[i]; 458 nz = ai[i+1] - ai[i]; 459 for (j=0; j<nz; j++) 460 rs += PetscAbsScalar(v[j]); 461 if (rs>sctx.shift_top) sctx.shift_top = rs; 462 } 463 sctx.shift_top *= 1.1; 464 sctx.nshift_max = 5; 465 sctx.shift_lo = 0.; 466 sctx.shift_hi = 1.; 467 } 468 469 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 470 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 471 ierr = PetscMalloc((n+1)*sizeof(MatScalar),&rtmp);CHKERRQ(ierr); 472 ics = ic; 473 474 do { 475 sctx.newshift = PETSC_FALSE; 476 for (i=0; i<n; i++) { 477 /* zero rtmp */ 478 /* L part */ 479 nz = bi[i+1] - bi[i]; 480 bjtmp = bj + bi[i]; 481 for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0; 482 483 /* U part */ 484 nz = bdiag[i]-bdiag[i+1]; 485 bjtmp = bj + bdiag[i+1]+1; 486 for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0; 487 488 /* load in initial (unfactored row) */ 489 nz = ai[r[i]+1] - ai[r[i]]; 490 ajtmp = aj + ai[r[i]]; 491 v = aa + ai[r[i]]; 492 for (j=0; j<nz; j++) { 493 rtmp[ics[ajtmp[j]]] = v[j]; 494 } 495 /* ZeropivotApply() */ 496 rtmp[i] += sctx.shift_amount; /* shift the diagonal of the matrix */ 497 498 /* elimination */ 499 bjtmp = bj + bi[i]; 500 row = *bjtmp++; 501 nzL = bi[i+1] - bi[i]; 502 for (k=0; k < nzL;k++) { 503 pc = rtmp + row; 504 if (*pc != (PetscScalar)0.0) { 505 pv = b->a + bdiag[row]; 506 multiplier = *pc * (*pv); 507 *pc = multiplier; 508 pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */ 509 pv = b->a + bdiag[row+1]+1; 510 nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */ 511 for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j]; 512 ierr = PetscLogFlops(2.0*nz);CHKERRQ(ierr); 513 } 514 row = *bjtmp++; 515 } 516 517 /* finished row so stick it into b->a */ 518 rs = 0.0; 519 /* L part */ 520 pv = b->a + bi[i] ; 521 pj = b->j + bi[i] ; 522 nz = bi[i+1] - bi[i]; 523 for (j=0; j<nz; j++) { 524 pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]); 525 } 526 527 /* U part */ 528 pv = b->a + bdiag[i+1]+1; 529 pj = b->j + bdiag[i+1]+1; 530 nz = bdiag[i] - bdiag[i+1]-1; 531 for (j=0; j<nz; j++) { 532 pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]); 533 } 534 535 sctx.rs = rs; 536 sctx.pv = rtmp[i]; 537 ierr = MatPivotCheck(A,info,&sctx,i);CHKERRQ(ierr); 538 if (sctx.newshift) break; /* break for-loop */ 539 rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */ 540 541 /* Mark diagonal and invert diagonal for simplier triangular solves */ 542 pv = b->a + bdiag[i]; 543 *pv = (PetscScalar)1.0/rtmp[i]; 544 545 } /* endof for (i=0; i<n; i++) { */ 546 547 /* MatPivotRefine() */ 548 if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) { 549 /* 550 * if no shift in this attempt & shifting & started shifting & can refine, 551 * then try lower shift 552 */ 553 sctx.shift_hi = sctx.shift_fraction; 554 sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.; 555 sctx.shift_amount = sctx.shift_fraction * sctx.shift_top; 556 sctx.newshift = PETSC_TRUE; 557 sctx.nshift++; 558 } 559 } while (sctx.newshift); 560 561 ierr = PetscFree(rtmp);CHKERRQ(ierr); 562 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 563 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 564 565 ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr); 566 ierr = ISIdentity(isicol,&col_identity);CHKERRQ(ierr); 567 if (row_identity && col_identity) { 568 C->ops->solve = MatSolve_SeqBAIJ_1_NaturalOrdering; 569 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering; 570 } else { 571 C->ops->solve = MatSolve_SeqBAIJ_1; 572 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1; 573 } 574 C->assembled = PETSC_TRUE; 575 ierr = PetscLogFlops(C->cmap->n);CHKERRQ(ierr); 576 577 /* MatShiftView(A,info,&sctx) */ 578 if (sctx.nshift) { 579 if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { 580 ierr = PetscInfo4(A,"number of shift_pd tries %D, shift_amount %G, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,sctx.shift_amount,sctx.shift_fraction,sctx.shift_top);CHKERRQ(ierr); 581 } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) { 582 ierr = PetscInfo2(A,"number of shift_nz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);CHKERRQ(ierr); 583 } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) { 584 ierr = PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %G\n",sctx.nshift,info->shiftamount);CHKERRQ(ierr); 585 } 586 } 587 PetscFunctionReturn(0); 588 } 589 590 #undef __FUNCT__ 591 #define __FUNCT__ "MatLUFactorNumeric_SeqBAIJ_1_inplace" 592 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo *info) 593 { 594 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ *)C->data; 595 IS isrow = b->row,isicol = b->icol; 596 PetscErrorCode ierr; 597 const PetscInt *r,*ic; 598 PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j; 599 PetscInt *ajtmpold,*ajtmp,nz,row,*ai = a->i,*aj = a->j; 600 PetscInt *diag_offset = b->diag,diag,*pj; 601 MatScalar *pv,*v,*rtmp,multiplier,*pc; 602 MatScalar *ba = b->a,*aa = a->a; 603 PetscBool row_identity, col_identity; 604 605 PetscFunctionBegin; 606 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 607 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 608 ierr = PetscMalloc((n+1)*sizeof(MatScalar),&rtmp);CHKERRQ(ierr); 609 610 for (i=0; i<n; i++) { 611 nz = bi[i+1] - bi[i]; 612 ajtmp = bj + bi[i]; 613 for (j=0; j<nz; j++) rtmp[ajtmp[j]] = 0.0; 614 615 /* load in initial (unfactored row) */ 616 nz = ai[r[i]+1] - ai[r[i]]; 617 ajtmpold = aj + ai[r[i]]; 618 v = aa + ai[r[i]]; 619 for (j=0; j<nz; j++) rtmp[ic[ajtmpold[j]]] = v[j]; 620 621 row = *ajtmp++; 622 while (row < i) { 623 pc = rtmp + row; 624 if (*pc != 0.0) { 625 pv = ba + diag_offset[row]; 626 pj = bj + diag_offset[row] + 1; 627 multiplier = *pc * *pv++; 628 *pc = multiplier; 629 nz = bi[row+1] - diag_offset[row] - 1; 630 for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j]; 631 ierr = PetscLogFlops(1.0+2.0*nz);CHKERRQ(ierr); 632 } 633 row = *ajtmp++; 634 } 635 /* finished row so stick it into b->a */ 636 pv = ba + bi[i]; 637 pj = bj + bi[i]; 638 nz = bi[i+1] - bi[i]; 639 for (j=0; j<nz; j++) {pv[j] = rtmp[pj[j]];} 640 diag = diag_offset[i] - bi[i]; 641 /* check pivot entry for current row */ 642 if (pv[diag] == 0.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot: row in original ordering %D in permuted ordering %D",r[i],i); 643 pv[diag] = 1.0/pv[diag]; 644 } 645 646 ierr = PetscFree(rtmp);CHKERRQ(ierr); 647 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 648 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 649 ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr); 650 ierr = ISIdentity(isicol,&col_identity);CHKERRQ(ierr); 651 if (row_identity && col_identity) { 652 C->ops->solve = MatSolve_SeqBAIJ_1_NaturalOrdering_inplace; 653 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace; 654 } else { 655 C->ops->solve = MatSolve_SeqBAIJ_1_inplace; 656 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_inplace; 657 } 658 C->assembled = PETSC_TRUE; 659 ierr = PetscLogFlops(C->cmap->n);CHKERRQ(ierr); 660 PetscFunctionReturn(0); 661 } 662 663 EXTERN_C_BEGIN 664 #undef __FUNCT__ 665 #define __FUNCT__ "MatGetFactor_seqbaij_petsc" 666 PetscErrorCode MatGetFactor_seqbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 667 { 668 PetscInt n = A->rmap->n; 669 PetscErrorCode ierr; 670 671 PetscFunctionBegin; 672 #if defined(PETSC_USE_COMPLEX) 673 if (A->hermitian && (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC))SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported"); 674 #endif 675 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 676 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 677 if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) { 678 ierr = MatSetType(*B,MATSEQBAIJ);CHKERRQ(ierr); 679 (*B)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqBAIJ; 680 (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqBAIJ; 681 } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 682 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 683 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 684 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqBAIJ; 685 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqBAIJ; 686 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 687 (*B)->factortype = ftype; 688 PetscFunctionReturn(0); 689 } 690 EXTERN_C_END 691 692 EXTERN_C_BEGIN 693 #undef __FUNCT__ 694 #define __FUNCT__ "MatGetFactorAvailable_seqbaij_petsc" 695 PetscErrorCode MatGetFactorAvailable_seqbaij_petsc(Mat A,MatFactorType ftype,PetscBool *flg) 696 { 697 PetscFunctionBegin; 698 *flg = PETSC_TRUE; 699 PetscFunctionReturn(0); 700 } 701 EXTERN_C_END 702 703 /* ----------------------------------------------------------- */ 704 #undef __FUNCT__ 705 #define __FUNCT__ "MatLUFactor_SeqBAIJ" 706 PetscErrorCode MatLUFactor_SeqBAIJ(Mat A,IS row,IS col,const MatFactorInfo *info) 707 { 708 PetscErrorCode ierr; 709 Mat C; 710 711 PetscFunctionBegin; 712 ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&C);CHKERRQ(ierr); 713 ierr = MatLUFactorSymbolic(C,A,row,col,info);CHKERRQ(ierr); 714 ierr = MatLUFactorNumeric(C,A,info);CHKERRQ(ierr); 715 A->ops->solve = C->ops->solve; 716 A->ops->solvetranspose = C->ops->solvetranspose; 717 ierr = MatHeaderMerge(A,C);CHKERRQ(ierr); 718 ierr = PetscLogObjectParent(A,((Mat_SeqBAIJ*)(A->data))->icol);CHKERRQ(ierr); 719 PetscFunctionReturn(0); 720 } 721 722 #include <../src/mat/impls/sbaij/seq/sbaij.h> 723 #undef __FUNCT__ 724 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqBAIJ_N" 725 PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N(Mat C,Mat A,const MatFactorInfo *info) 726 { 727 PetscErrorCode ierr; 728 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 729 Mat_SeqSBAIJ *b=(Mat_SeqSBAIJ*)C->data; 730 IS ip=b->row; 731 const PetscInt *rip; 732 PetscInt i,j,mbs=a->mbs,bs=A->rmap->bs,*bi=b->i,*bj=b->j,*bcol; 733 PetscInt *ai=a->i,*aj=a->j; 734 PetscInt k,jmin,jmax,*jl,*il,col,nexti,ili,nz; 735 MatScalar *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi; 736 PetscReal rs; 737 FactorShiftCtx sctx; 738 739 PetscFunctionBegin; 740 if (bs > 1) {/* convert A to a SBAIJ matrix and apply Cholesky factorization from it */ 741 if (!a->sbaijMat) { 742 ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr); 743 } 744 ierr = (a->sbaijMat)->ops->choleskyfactornumeric(C,a->sbaijMat,info);CHKERRQ(ierr); 745 ierr = MatDestroy(&a->sbaijMat);CHKERRQ(ierr); 746 PetscFunctionReturn(0); 747 } 748 749 /* MatPivotSetUp(): initialize shift context sctx */ 750 ierr = PetscMemzero(&sctx,sizeof(FactorShiftCtx));CHKERRQ(ierr); 751 752 ierr = ISGetIndices(ip,&rip);CHKERRQ(ierr); 753 ierr = PetscMalloc3(mbs,MatScalar,&rtmp,mbs,PetscInt,&il,mbs,PetscInt,&jl);CHKERRQ(ierr); 754 755 sctx.shift_amount = 0.; 756 sctx.nshift = 0; 757 do { 758 sctx.newshift = PETSC_FALSE; 759 for (i=0; i<mbs; i++) { 760 rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0; 761 } 762 763 for (k = 0; k<mbs; k++) { 764 bval = ba + bi[k]; 765 /* initialize k-th row by the perm[k]-th row of A */ 766 jmin = ai[rip[k]]; jmax = ai[rip[k]+1]; 767 for (j = jmin; j < jmax; j++) { 768 col = rip[aj[j]]; 769 if (col >= k) { /* only take upper triangular entry */ 770 rtmp[col] = aa[j]; 771 *bval++ = 0.0; /* for in-place factorization */ 772 } 773 } 774 775 /* shift the diagonal of the matrix */ 776 if (sctx.nshift) rtmp[k] += sctx.shift_amount; 777 778 /* modify k-th row by adding in those rows i with U(i,k)!=0 */ 779 dk = rtmp[k]; 780 i = jl[k]; /* first row to be added to k_th row */ 781 782 while (i < k) { 783 nexti = jl[i]; /* next row to be added to k_th row */ 784 785 /* compute multiplier, update diag(k) and U(i,k) */ 786 ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */ 787 uikdi = - ba[ili]*ba[bi[i]]; /* diagonal(k) */ 788 dk += uikdi*ba[ili]; 789 ba[ili] = uikdi; /* -U(i,k) */ 790 791 /* add multiple of row i to k-th row */ 792 jmin = ili + 1; jmax = bi[i+1]; 793 if (jmin < jmax) { 794 for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j]; 795 /* update il and jl for row i */ 796 il[i] = jmin; 797 j = bj[jmin]; jl[i] = jl[j]; jl[j] = i; 798 } 799 i = nexti; 800 } 801 802 /* shift the diagonals when zero pivot is detected */ 803 /* compute rs=sum of abs(off-diagonal) */ 804 rs = 0.0; 805 jmin = bi[k]+1; 806 nz = bi[k+1] - jmin; 807 if (nz) { 808 bcol = bj + jmin; 809 while (nz--) { 810 rs += PetscAbsScalar(rtmp[*bcol]); 811 bcol++; 812 } 813 } 814 815 sctx.rs = rs; 816 sctx.pv = dk; 817 ierr = MatPivotCheck(A,info,&sctx,k);CHKERRQ(ierr); 818 if (sctx.newshift) break; 819 dk = sctx.pv; 820 821 /* copy data into U(k,:) */ 822 ba[bi[k]] = 1.0/dk; /* U(k,k) */ 823 jmin = bi[k]+1; jmax = bi[k+1]; 824 if (jmin < jmax) { 825 for (j=jmin; j<jmax; j++) { 826 col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0; 827 } 828 /* add the k-th row into il and jl */ 829 il[k] = jmin; 830 i = bj[jmin]; jl[k] = jl[i]; jl[i] = k; 831 } 832 } 833 } while (sctx.newshift); 834 ierr = PetscFree3(rtmp,il,jl);CHKERRQ(ierr); 835 836 ierr = ISRestoreIndices(ip,&rip);CHKERRQ(ierr); 837 C->assembled = PETSC_TRUE; 838 C->preallocated = PETSC_TRUE; 839 ierr = PetscLogFlops(C->rmap->N);CHKERRQ(ierr); 840 if (sctx.nshift) { 841 if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { 842 ierr = PetscInfo2(A,"number of shiftpd tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);CHKERRQ(ierr); 843 } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) { 844 ierr = PetscInfo2(A,"number of shiftnz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);CHKERRQ(ierr); 845 } 846 } 847 PetscFunctionReturn(0); 848 } 849 850 #undef __FUNCT__ 851 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering" 852 PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info) 853 { 854 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 855 Mat_SeqSBAIJ *b=(Mat_SeqSBAIJ*)C->data; 856 PetscErrorCode ierr; 857 PetscInt i,j,am=a->mbs; 858 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j; 859 PetscInt k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz; 860 MatScalar *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval; 861 PetscReal rs; 862 FactorShiftCtx sctx; 863 864 PetscFunctionBegin; 865 /* MatPivotSetUp(): initialize shift context sctx */ 866 ierr = PetscMemzero(&sctx,sizeof(FactorShiftCtx));CHKERRQ(ierr); 867 868 ierr = PetscMalloc3(am,MatScalar,&rtmp,am,PetscInt,&il,am,PetscInt,&jl);CHKERRQ(ierr); 869 870 do { 871 sctx.newshift = PETSC_FALSE; 872 for (i=0; i<am; i++) { 873 rtmp[i] = 0.0; jl[i] = am; il[0] = 0; 874 } 875 876 for (k = 0; k<am; k++) { 877 /* initialize k-th row with elements nonzero in row perm(k) of A */ 878 nz = ai[k+1] - ai[k]; 879 acol = aj + ai[k]; 880 aval = aa + ai[k]; 881 bval = ba + bi[k]; 882 while (nz--) { 883 if (*acol < k) { /* skip lower triangular entries */ 884 acol++; aval++; 885 } else { 886 rtmp[*acol++] = *aval++; 887 *bval++ = 0.0; /* for in-place factorization */ 888 } 889 } 890 891 /* shift the diagonal of the matrix */ 892 if (sctx.nshift) rtmp[k] += sctx.shift_amount; 893 894 /* modify k-th row by adding in those rows i with U(i,k)!=0 */ 895 dk = rtmp[k]; 896 i = jl[k]; /* first row to be added to k_th row */ 897 898 while (i < k) { 899 nexti = jl[i]; /* next row to be added to k_th row */ 900 /* compute multiplier, update D(k) and U(i,k) */ 901 ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */ 902 uikdi = - ba[ili]*ba[bi[i]]; 903 dk += uikdi*ba[ili]; 904 ba[ili] = uikdi; /* -U(i,k) */ 905 906 /* add multiple of row i to k-th row ... */ 907 jmin = ili + 1; 908 nz = bi[i+1] - jmin; 909 if (nz > 0) { 910 bcol = bj + jmin; 911 bval = ba + jmin; 912 while (nz --) rtmp[*bcol++] += uikdi*(*bval++); 913 /* update il and jl for i-th row */ 914 il[i] = jmin; 915 j = bj[jmin]; jl[i] = jl[j]; jl[j] = i; 916 } 917 i = nexti; 918 } 919 920 /* shift the diagonals when zero pivot is detected */ 921 /* compute rs=sum of abs(off-diagonal) */ 922 rs = 0.0; 923 jmin = bi[k]+1; 924 nz = bi[k+1] - jmin; 925 if (nz) { 926 bcol = bj + jmin; 927 while (nz--) { 928 rs += PetscAbsScalar(rtmp[*bcol]); 929 bcol++; 930 } 931 } 932 933 sctx.rs = rs; 934 sctx.pv = dk; 935 ierr = MatPivotCheck(A,info,&sctx,k);CHKERRQ(ierr); 936 if (sctx.newshift) break; /* sctx.shift_amount is updated */ 937 dk = sctx.pv; 938 939 /* copy data into U(k,:) */ 940 ba[bi[k]] = 1.0/dk; 941 jmin = bi[k]+1; 942 nz = bi[k+1] - jmin; 943 if (nz) { 944 bcol = bj + jmin; 945 bval = ba + jmin; 946 while (nz--) { 947 *bval++ = rtmp[*bcol]; 948 rtmp[*bcol++] = 0.0; 949 } 950 /* add k-th row into il and jl */ 951 il[k] = jmin; 952 i = bj[jmin]; jl[k] = jl[i]; jl[i] = k; 953 } 954 } 955 } while (sctx.newshift); 956 ierr = PetscFree3(rtmp,il,jl);CHKERRQ(ierr); 957 958 C->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 959 C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 960 C->assembled = PETSC_TRUE; 961 C->preallocated = PETSC_TRUE; 962 ierr = PetscLogFlops(C->rmap->N);CHKERRQ(ierr); 963 if (sctx.nshift) { 964 if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) { 965 ierr = PetscInfo2(A,"number of shiftnz tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);CHKERRQ(ierr); 966 } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { 967 ierr = PetscInfo2(A,"number of shiftpd tries %D, shift_amount %G\n",sctx.nshift,sctx.shift_amount);CHKERRQ(ierr); 968 } 969 } 970 PetscFunctionReturn(0); 971 } 972 973 #include <petscbt.h> 974 #include <../src/mat/utils/freespace.h> 975 #undef __FUNCT__ 976 #define __FUNCT__ "MatICCFactorSymbolic_SeqBAIJ" 977 PetscErrorCode MatICCFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info) 978 { 979 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 980 Mat_SeqSBAIJ *b; 981 Mat B; 982 PetscErrorCode ierr; 983 PetscBool perm_identity; 984 PetscInt reallocs=0,i,*ai=a->i,*aj=a->j,am=a->mbs,bs=A->rmap->bs,*ui; 985 const PetscInt *rip; 986 PetscInt jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow; 987 PetscInt nlnk,*lnk,*lnk_lvl=PETSC_NULL,ncols,ncols_upper,*cols,*cols_lvl,*uj,**uj_ptr,**uj_lvl_ptr; 988 PetscReal fill=info->fill,levels=info->levels; 989 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 990 PetscFreeSpaceList free_space_lvl=PETSC_NULL,current_space_lvl=PETSC_NULL; 991 PetscBT lnkbt; 992 993 PetscFunctionBegin; 994 if (bs > 1) { 995 if (!a->sbaijMat) { 996 ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr); 997 } 998 (fact)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */ 999 ierr = MatICCFactorSymbolic(fact,a->sbaijMat,perm,info);CHKERRQ(ierr); 1000 PetscFunctionReturn(0); 1001 } 1002 1003 ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr); 1004 ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr); 1005 1006 /* special case that simply copies fill pattern */ 1007 if (!levels && perm_identity) { 1008 ierr = MatMarkDiagonal_SeqBAIJ(A);CHKERRQ(ierr); 1009 ierr = PetscMalloc((am+1)*sizeof(PetscInt),&ui);CHKERRQ(ierr); 1010 for (i=0; i<am; i++) { 1011 ui[i] = ai[i+1] - a->diag[i]; /* ui: rowlengths - changes when !perm_identity */ 1012 } 1013 B = fact; 1014 ierr = MatSeqSBAIJSetPreallocation(B,1,0,ui);CHKERRQ(ierr); 1015 1016 1017 b = (Mat_SeqSBAIJ*)B->data; 1018 uj = b->j; 1019 for (i=0; i<am; i++) { 1020 aj = a->j + a->diag[i]; 1021 for (j=0; j<ui[i]; j++) { 1022 *uj++ = *aj++; 1023 } 1024 b->ilen[i] = ui[i]; 1025 } 1026 ierr = PetscFree(ui);CHKERRQ(ierr); 1027 B->factortype = MAT_FACTOR_NONE; 1028 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1029 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1030 B->factortype = MAT_FACTOR_ICC; 1031 1032 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering; 1033 PetscFunctionReturn(0); 1034 } 1035 1036 /* initialization */ 1037 ierr = PetscMalloc((am+1)*sizeof(PetscInt),&ui);CHKERRQ(ierr); 1038 ui[0] = 0; 1039 ierr = PetscMalloc((2*am+1)*sizeof(PetscInt),&cols_lvl);CHKERRQ(ierr); 1040 1041 /* jl: linked list for storing indices of the pivot rows 1042 il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */ 1043 ierr = PetscMalloc4(am,PetscInt*,&uj_ptr,am,PetscInt*,&uj_lvl_ptr,am,PetscInt,&il,am,PetscInt,&jl);CHKERRQ(ierr); 1044 for (i=0; i<am; i++) { 1045 jl[i] = am; il[i] = 0; 1046 } 1047 1048 /* create and initialize a linked list for storing column indices of the active row k */ 1049 nlnk = am + 1; 1050 ierr = PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr); 1051 1052 /* initial FreeSpace size is fill*(ai[am]+am)/2 */ 1053 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+am)/2),&free_space);CHKERRQ(ierr); 1054 current_space = free_space; 1055 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+am)/2),&free_space_lvl);CHKERRQ(ierr); 1056 current_space_lvl = free_space_lvl; 1057 1058 for (k=0; k<am; k++) { /* for each active row k */ 1059 /* initialize lnk by the column indices of row rip[k] of A */ 1060 nzk = 0; 1061 ncols = ai[rip[k]+1] - ai[rip[k]]; 1062 ncols_upper = 0; 1063 cols = cols_lvl + am; 1064 for (j=0; j<ncols; j++) { 1065 i = rip[*(aj + ai[rip[k]] + j)]; 1066 if (i >= k) { /* only take upper triangular entry */ 1067 cols[ncols_upper] = i; 1068 cols_lvl[ncols_upper] = -1; /* initialize level for nonzero entries */ 1069 ncols_upper++; 1070 } 1071 } 1072 ierr = PetscIncompleteLLAdd(ncols_upper,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr); 1073 nzk += nlnk; 1074 1075 /* update lnk by computing fill-in for each pivot row to be merged in */ 1076 prow = jl[k]; /* 1st pivot row */ 1077 1078 while (prow < k) { 1079 nextprow = jl[prow]; 1080 1081 /* merge prow into k-th row */ 1082 jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */ 1083 jmax = ui[prow+1]; 1084 ncols = jmax-jmin; 1085 i = jmin - ui[prow]; 1086 cols = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */ 1087 for (j=0; j<ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j); 1088 ierr = PetscIncompleteLLAddSorted(ncols,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr); 1089 nzk += nlnk; 1090 1091 /* update il and jl for prow */ 1092 if (jmin < jmax) { 1093 il[prow] = jmin; 1094 j = *cols; jl[prow] = jl[j]; jl[j] = prow; 1095 } 1096 prow = nextprow; 1097 } 1098 1099 /* if free space is not available, make more free space */ 1100 if (current_space->local_remaining<nzk) { 1101 i = am - k + 1; /* num of unfactored rows */ 1102 i = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */ 1103 ierr = PetscFreeSpaceGet(i,¤t_space);CHKERRQ(ierr); 1104 ierr = PetscFreeSpaceGet(i,¤t_space_lvl);CHKERRQ(ierr); 1105 reallocs++; 1106 } 1107 1108 /* copy data into free_space and free_space_lvl, then initialize lnk */ 1109 ierr = PetscIncompleteLLClean(am,am,nzk,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);CHKERRQ(ierr); 1110 1111 /* add the k-th row into il and jl */ 1112 if (nzk-1 > 0) { 1113 i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */ 1114 jl[k] = jl[i]; jl[i] = k; 1115 il[k] = ui[k] + 1; 1116 } 1117 uj_ptr[k] = current_space->array; 1118 uj_lvl_ptr[k] = current_space_lvl->array; 1119 1120 current_space->array += nzk; 1121 current_space->local_used += nzk; 1122 current_space->local_remaining -= nzk; 1123 1124 current_space_lvl->array += nzk; 1125 current_space_lvl->local_used += nzk; 1126 current_space_lvl->local_remaining -= nzk; 1127 1128 ui[k+1] = ui[k] + nzk; 1129 } 1130 1131 ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr); 1132 ierr = PetscFree4(uj_ptr,uj_lvl_ptr,il,jl);CHKERRQ(ierr); 1133 ierr = PetscFree(cols_lvl);CHKERRQ(ierr); 1134 1135 /* copy free_space into uj and free free_space; set uj in new datastructure; */ 1136 ierr = PetscMalloc((ui[am]+1)*sizeof(PetscInt),&uj);CHKERRQ(ierr); 1137 ierr = PetscFreeSpaceContiguous(&free_space,uj);CHKERRQ(ierr); 1138 ierr = PetscIncompleteLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 1139 ierr = PetscFreeSpaceDestroy(free_space_lvl);CHKERRQ(ierr); 1140 1141 /* put together the new matrix in MATSEQSBAIJ format */ 1142 B = fact; 1143 ierr = MatSeqSBAIJSetPreallocation(B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1144 1145 b = (Mat_SeqSBAIJ*)B->data; 1146 b->singlemalloc = PETSC_FALSE; 1147 b->free_a = PETSC_TRUE; 1148 b->free_ij = PETSC_TRUE; 1149 ierr = PetscMalloc((ui[am]+1)*sizeof(MatScalar),&b->a);CHKERRQ(ierr); 1150 b->j = uj; 1151 b->i = ui; 1152 b->diag = 0; 1153 b->ilen = 0; 1154 b->imax = 0; 1155 b->row = perm; 1156 b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */ 1157 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 1158 b->icol = perm; 1159 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 1160 ierr = PetscMalloc((am+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr); 1161 ierr = PetscLogObjectMemory(B,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr); 1162 b->maxnz = b->nz = ui[am]; 1163 1164 B->info.factor_mallocs = reallocs; 1165 B->info.fill_ratio_given = fill; 1166 if (ai[am] != 0.) { 1167 /* nonzeros in lower triangular part of A (includign diagonals)= (ai[am]+am)/2 */ 1168 B->info.fill_ratio_needed = ((PetscReal)2*ui[am])/(ai[am]+am); 1169 } else { 1170 B->info.fill_ratio_needed = 0.0; 1171 } 1172 #if defined(PETSC_USE_INFO) 1173 if (ai[am] != 0) { 1174 PetscReal af = B->info.fill_ratio_needed; 1175 ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,fill,af);CHKERRQ(ierr); 1176 ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr); 1177 ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G) for best performance.\n",af);CHKERRQ(ierr); 1178 } else { 1179 ierr = PetscInfo(A,"Empty matrix.\n");CHKERRQ(ierr); 1180 } 1181 #endif 1182 if (perm_identity) { 1183 B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 1184 B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 1185 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering; 1186 } else { 1187 (fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N; 1188 } 1189 PetscFunctionReturn(0); 1190 } 1191 1192 #undef __FUNCT__ 1193 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqBAIJ" 1194 PetscErrorCode MatCholeskyFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info) 1195 { 1196 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1197 Mat_SeqSBAIJ *b; 1198 Mat B; 1199 PetscErrorCode ierr; 1200 PetscBool perm_identity; 1201 PetscReal fill = info->fill; 1202 const PetscInt *rip; 1203 PetscInt i,mbs=a->mbs,bs=A->rmap->bs,*ai=a->i,*aj=a->j,reallocs=0,prow; 1204 PetscInt *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow; 1205 PetscInt nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr; 1206 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 1207 PetscBT lnkbt; 1208 1209 PetscFunctionBegin; 1210 if (bs > 1) { /* convert to seqsbaij */ 1211 if (!a->sbaijMat) { 1212 ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr); 1213 } 1214 (fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */ 1215 ierr = MatCholeskyFactorSymbolic(fact,a->sbaijMat,perm,info);CHKERRQ(ierr); 1216 PetscFunctionReturn(0); 1217 } 1218 1219 /* check whether perm is the identity mapping */ 1220 ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr); 1221 if (!perm_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1222 ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr); 1223 1224 /* initialization */ 1225 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&ui);CHKERRQ(ierr); 1226 ui[0] = 0; 1227 1228 /* jl: linked list for storing indices of the pivot rows 1229 il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */ 1230 ierr = PetscMalloc4(mbs,PetscInt*,&ui_ptr,mbs,PetscInt,&il,mbs,PetscInt,&jl,mbs,PetscInt,&cols);CHKERRQ(ierr); 1231 for (i=0; i<mbs; i++) { 1232 jl[i] = mbs; il[i] = 0; 1233 } 1234 1235 /* create and initialize a linked list for storing column indices of the active row k */ 1236 nlnk = mbs + 1; 1237 ierr = PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1238 1239 /* initial FreeSpace size is fill* (ai[mbs]+mbs)/2 */ 1240 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[mbs]+mbs)/2),&free_space);CHKERRQ(ierr); 1241 current_space = free_space; 1242 1243 for (k=0; k<mbs; k++) { /* for each active row k */ 1244 /* initialize lnk by the column indices of row rip[k] of A */ 1245 nzk = 0; 1246 ncols = ai[rip[k]+1] - ai[rip[k]]; 1247 ncols_upper = 0; 1248 for (j=0; j<ncols; j++) { 1249 i = rip[*(aj + ai[rip[k]] + j)]; 1250 if (i >= k) { /* only take upper triangular entry */ 1251 cols[ncols_upper] = i; 1252 ncols_upper++; 1253 } 1254 } 1255 ierr = PetscLLAdd(ncols_upper,cols,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1256 nzk += nlnk; 1257 1258 /* update lnk by computing fill-in for each pivot row to be merged in */ 1259 prow = jl[k]; /* 1st pivot row */ 1260 1261 while (prow < k) { 1262 nextprow = jl[prow]; 1263 /* merge prow into k-th row */ 1264 jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */ 1265 jmax = ui[prow+1]; 1266 ncols = jmax-jmin; 1267 uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */ 1268 ierr = PetscLLAddSorted(ncols,uj_ptr,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1269 nzk += nlnk; 1270 1271 /* update il and jl for prow */ 1272 if (jmin < jmax) { 1273 il[prow] = jmin; 1274 j = *uj_ptr; jl[prow] = jl[j]; jl[j] = prow; 1275 } 1276 prow = nextprow; 1277 } 1278 1279 /* if free space is not available, make more free space */ 1280 if (current_space->local_remaining<nzk) { 1281 i = mbs - k + 1; /* num of unfactored rows */ 1282 i = PetscMin(i*nzk, i*(i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */ 1283 ierr = PetscFreeSpaceGet(i,¤t_space);CHKERRQ(ierr); 1284 reallocs++; 1285 } 1286 1287 /* copy data into free space, then initialize lnk */ 1288 ierr = PetscLLClean(mbs,mbs,nzk,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 1289 1290 /* add the k-th row into il and jl */ 1291 if (nzk-1 > 0) { 1292 i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */ 1293 jl[k] = jl[i]; jl[i] = k; 1294 il[k] = ui[k] + 1; 1295 } 1296 ui_ptr[k] = current_space->array; 1297 current_space->array += nzk; 1298 current_space->local_used += nzk; 1299 current_space->local_remaining -= nzk; 1300 1301 ui[k+1] = ui[k] + nzk; 1302 } 1303 1304 ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr); 1305 ierr = PetscFree4(ui_ptr,il,jl,cols);CHKERRQ(ierr); 1306 1307 /* copy free_space into uj and free free_space; set uj in new datastructure; */ 1308 ierr = PetscMalloc((ui[mbs]+1)*sizeof(PetscInt),&uj);CHKERRQ(ierr); 1309 ierr = PetscFreeSpaceContiguous(&free_space,uj);CHKERRQ(ierr); 1310 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 1311 1312 /* put together the new matrix in MATSEQSBAIJ format */ 1313 B = fact; 1314 ierr = MatSeqSBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1315 1316 b = (Mat_SeqSBAIJ*)B->data; 1317 b->singlemalloc = PETSC_FALSE; 1318 b->free_a = PETSC_TRUE; 1319 b->free_ij = PETSC_TRUE; 1320 ierr = PetscMalloc((ui[mbs]+1)*sizeof(MatScalar),&b->a);CHKERRQ(ierr); 1321 b->j = uj; 1322 b->i = ui; 1323 b->diag = 0; 1324 b->ilen = 0; 1325 b->imax = 0; 1326 b->row = perm; 1327 b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */ 1328 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 1329 b->icol = perm; 1330 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 1331 ierr = PetscMalloc((mbs+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr); 1332 ierr = PetscLogObjectMemory(B,(ui[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr); 1333 b->maxnz = b->nz = ui[mbs]; 1334 1335 B->info.factor_mallocs = reallocs; 1336 B->info.fill_ratio_given = fill; 1337 if (ai[mbs] != 0.) { 1338 /* nonzeros in lower triangular part of A = (ai[mbs]+mbs)/2 */ 1339 B->info.fill_ratio_needed = ((PetscReal)2*ui[mbs])/(ai[mbs]+mbs); 1340 } else { 1341 B->info.fill_ratio_needed = 0.0; 1342 } 1343 #if defined(PETSC_USE_INFO) 1344 if (ai[mbs] != 0.) { 1345 PetscReal af = B->info.fill_ratio_needed; 1346 ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,fill,af);CHKERRQ(ierr); 1347 ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr); 1348 ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G) for best performance.\n",af);CHKERRQ(ierr); 1349 } else { 1350 ierr = PetscInfo(A,"Empty matrix.\n");CHKERRQ(ierr); 1351 } 1352 #endif 1353 if (perm_identity) { 1354 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering; 1355 } else { 1356 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N; 1357 } 1358 PetscFunctionReturn(0); 1359 } 1360 1361 #undef __FUNCT__ 1362 #define __FUNCT__ "MatSolve_SeqBAIJ_N_NaturalOrdering" 1363 PetscErrorCode MatSolve_SeqBAIJ_N_NaturalOrdering(Mat A,Vec bb,Vec xx) 1364 { 1365 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 1366 PetscErrorCode ierr; 1367 const PetscInt *ai=a->i,*aj=a->j,*adiag=a->diag,*vi; 1368 PetscInt i,k,n=a->mbs; 1369 PetscInt nz,bs=A->rmap->bs,bs2=a->bs2; 1370 MatScalar *aa=a->a,*v; 1371 PetscScalar *x,*b,*s,*t,*ls; 1372 1373 PetscFunctionBegin; 1374 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 1375 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1376 t = a->solve_work; 1377 1378 /* forward solve the lower triangular */ 1379 ierr = PetscMemcpy(t,b,bs*sizeof(PetscScalar));CHKERRQ(ierr); /* copy 1st block of b to t */ 1380 1381 for (i=1; i<n; i++) { 1382 v = aa + bs2*ai[i]; 1383 vi = aj + ai[i]; 1384 nz = ai[i+1] - ai[i]; 1385 s = t + bs*i; 1386 ierr = PetscMemcpy(s,b+bs*i,bs*sizeof(PetscScalar));CHKERRQ(ierr); /* copy i_th block of b to t */ 1387 for (k=0;k<nz;k++) { 1388 PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[k]); 1389 v += bs2; 1390 } 1391 } 1392 1393 /* backward solve the upper triangular */ 1394 ls = a->solve_work + A->cmap->n; 1395 for (i=n-1; i>=0; i--) { 1396 v = aa + bs2*(adiag[i+1]+1); 1397 vi = aj + adiag[i+1]+1; 1398 nz = adiag[i] - adiag[i+1]-1; 1399 ierr = PetscMemcpy(ls,t+i*bs,bs*sizeof(PetscScalar));CHKERRQ(ierr); 1400 for (k=0;k<nz;k++) { 1401 PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[k]); 1402 v += bs2; 1403 } 1404 PetscKernel_w_gets_A_times_v(bs,ls,aa+bs2*adiag[i],t+i*bs); /* *inv(diagonal[i]) */ 1405 ierr = PetscMemcpy(x+i*bs,t+i*bs,bs*sizeof(PetscScalar));CHKERRQ(ierr); 1406 } 1407 1408 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 1409 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1410 ierr = PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);CHKERRQ(ierr); 1411 PetscFunctionReturn(0); 1412 } 1413 1414 #undef __FUNCT__ 1415 #define __FUNCT__ "MatSolve_SeqBAIJ_N" 1416 PetscErrorCode MatSolve_SeqBAIJ_N(Mat A,Vec bb,Vec xx) 1417 { 1418 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 1419 IS iscol=a->col,isrow=a->row; 1420 PetscErrorCode ierr; 1421 const PetscInt *r,*c,*rout,*cout,*ai=a->i,*aj=a->j,*adiag=a->diag,*vi; 1422 PetscInt i,m,n=a->mbs; 1423 PetscInt nz,bs=A->rmap->bs,bs2=a->bs2; 1424 MatScalar *aa=a->a,*v; 1425 PetscScalar *x,*b,*s,*t,*ls; 1426 1427 PetscFunctionBegin; 1428 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 1429 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1430 t = a->solve_work; 1431 1432 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 1433 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout; 1434 1435 /* forward solve the lower triangular */ 1436 ierr = PetscMemcpy(t,b+bs*r[0],bs*sizeof(PetscScalar));CHKERRQ(ierr); 1437 for (i=1; i<n; i++) { 1438 v = aa + bs2*ai[i]; 1439 vi = aj + ai[i]; 1440 nz = ai[i+1] - ai[i]; 1441 s = t + bs*i; 1442 ierr = PetscMemcpy(s,b+bs*r[i],bs*sizeof(PetscScalar));CHKERRQ(ierr); 1443 for (m=0;m<nz;m++) { 1444 PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[m]); 1445 v += bs2; 1446 } 1447 } 1448 1449 /* backward solve the upper triangular */ 1450 ls = a->solve_work + A->cmap->n; 1451 for (i=n-1; i>=0; i--) { 1452 v = aa + bs2*(adiag[i+1]+1); 1453 vi = aj + adiag[i+1]+1; 1454 nz = adiag[i] - adiag[i+1] - 1; 1455 ierr = PetscMemcpy(ls,t+i*bs,bs*sizeof(PetscScalar));CHKERRQ(ierr); 1456 for (m=0;m<nz;m++) { 1457 PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[m]); 1458 v += bs2; 1459 } 1460 PetscKernel_w_gets_A_times_v(bs,ls,v,t+i*bs); /* *inv(diagonal[i]) */ 1461 ierr = PetscMemcpy(x + bs*c[i],t+i*bs,bs*sizeof(PetscScalar));CHKERRQ(ierr); 1462 } 1463 ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr); 1464 ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr); 1465 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 1466 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1467 ierr = PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);CHKERRQ(ierr); 1468 PetscFunctionReturn(0); 1469 } 1470 1471 #undef __FUNCT__ 1472 #define __FUNCT__ "MatBlockAbs_privat" 1473 /* 1474 For each block in an block array saves the largest absolute value in the block into another array 1475 */ 1476 static PetscErrorCode MatBlockAbs_private(PetscInt nbs,PetscInt bs2,PetscScalar *blockarray,PetscReal *absarray) 1477 { 1478 PetscErrorCode ierr; 1479 PetscInt i,j; 1480 1481 PetscFunctionBegin; 1482 ierr = PetscMemzero(absarray,(nbs+1)*sizeof(PetscReal));CHKERRQ(ierr); 1483 for (i=0; i<nbs; i++) { 1484 for (j=0; j<bs2; j++) { 1485 if (absarray[i] < PetscAbsScalar(blockarray[i*nbs+j])) absarray[i] = PetscAbsScalar(blockarray[i*nbs+j]); 1486 } 1487 } 1488 PetscFunctionReturn(0); 1489 } 1490 1491 #undef __FUNCT__ 1492 #define __FUNCT__ "MatILUDTFactor_SeqBAIJ" 1493 /* 1494 This needs to be renamed and called by the regular MatILUFactor_SeqBAIJ when drop tolerance is used 1495 */ 1496 PetscErrorCode MatILUDTFactor_SeqBAIJ(Mat A,IS isrow,IS iscol,const MatFactorInfo *info,Mat *fact) 1497 { 1498 Mat B = *fact; 1499 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b; 1500 IS isicol; 1501 PetscErrorCode ierr; 1502 const PetscInt *r,*ic; 1503 PetscInt i,mbs=a->mbs,bs=A->rmap->bs,bs2=a->bs2,*ai=a->i,*aj=a->j,*ajtmp,*adiag; 1504 PetscInt *bi,*bj,*bdiag; 1505 1506 PetscInt row,nzi,nzi_bl,nzi_bu,*im,dtcount,nzi_al,nzi_au; 1507 PetscInt nlnk,*lnk; 1508 PetscBT lnkbt; 1509 PetscBool row_identity,icol_identity; 1510 MatScalar *aatmp,*pv,*batmp,*ba,*rtmp,*pc,*multiplier,*vtmp; 1511 PetscInt j,nz,*pj,*bjtmp,k,ncut,*jtmp; 1512 1513 PetscReal dt=info->dt; /* shift=info->shiftamount; */ 1514 PetscInt nnz_max; 1515 PetscBool missing; 1516 PetscReal *vtmp_abs; 1517 MatScalar *v_work; 1518 PetscInt *v_pivots; 1519 1520 PetscFunctionBegin; 1521 /* ------- symbolic factorization, can be reused ---------*/ 1522 ierr = MatMissingDiagonal(A,&missing,&i);CHKERRQ(ierr); 1523 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i); 1524 adiag=a->diag; 1525 1526 ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr); 1527 1528 /* bdiag is location of diagonal in factor */ 1529 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&bdiag);CHKERRQ(ierr); 1530 1531 /* allocate row pointers bi */ 1532 ierr = PetscMalloc((2*mbs+2)*sizeof(PetscInt),&bi);CHKERRQ(ierr); 1533 1534 /* allocate bj and ba; max num of nonzero entries is (ai[n]+2*n*dtcount+2) */ 1535 dtcount = (PetscInt)info->dtcount; 1536 if (dtcount > mbs-1) dtcount = mbs-1; 1537 nnz_max = ai[mbs]+2*mbs*dtcount +2; 1538 /* printf("MatILUDTFactor_SeqBAIJ, bs %d, ai[mbs] %d, nnz_max %d, dtcount %d\n",bs,ai[mbs],nnz_max,dtcount); */ 1539 ierr = PetscMalloc(nnz_max*sizeof(PetscInt),&bj);CHKERRQ(ierr); 1540 nnz_max = nnz_max*bs2; 1541 ierr = PetscMalloc(nnz_max*sizeof(MatScalar),&ba);CHKERRQ(ierr); 1542 1543 /* put together the new matrix */ 1544 ierr = MatSeqBAIJSetPreallocation_SeqBAIJ(B,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1545 ierr = PetscLogObjectParent(B,isicol);CHKERRQ(ierr); 1546 b = (Mat_SeqBAIJ*)(B)->data; 1547 b->free_a = PETSC_TRUE; 1548 b->free_ij = PETSC_TRUE; 1549 b->singlemalloc = PETSC_FALSE; 1550 b->a = ba; 1551 b->j = bj; 1552 b->i = bi; 1553 b->diag = bdiag; 1554 b->ilen = 0; 1555 b->imax = 0; 1556 b->row = isrow; 1557 b->col = iscol; 1558 ierr = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr); 1559 ierr = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr); 1560 b->icol = isicol; 1561 ierr = PetscMalloc((bs*(mbs+1))*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr); 1562 1563 ierr = PetscLogObjectMemory(B,nnz_max*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr); 1564 b->maxnz = nnz_max/bs2; 1565 1566 (B)->factortype = MAT_FACTOR_ILUDT; 1567 (B)->info.factor_mallocs = 0; 1568 (B)->info.fill_ratio_given = ((PetscReal)nnz_max)/((PetscReal)(ai[mbs]*bs2)); 1569 CHKMEMQ; 1570 /* ------- end of symbolic factorization ---------*/ 1571 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 1572 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 1573 1574 /* linked list for storing column indices of the active row */ 1575 nlnk = mbs + 1; 1576 ierr = PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1577 1578 /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */ 1579 ierr = PetscMalloc2(mbs,PetscInt,&im,mbs,PetscInt,&jtmp);CHKERRQ(ierr); 1580 /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */ 1581 ierr = PetscMalloc2(mbs*bs2,MatScalar,&rtmp,mbs*bs2,MatScalar,&vtmp);CHKERRQ(ierr); 1582 ierr = PetscMalloc((mbs+1)*sizeof(PetscReal),&vtmp_abs);CHKERRQ(ierr); 1583 ierr = PetscMalloc3(bs,MatScalar,&v_work,bs2,MatScalar,&multiplier,bs,PetscInt,&v_pivots);CHKERRQ(ierr); 1584 1585 bi[0] = 0; 1586 bdiag[0] = (nnz_max/bs2)-1; /* location of diagonal in factor B */ 1587 bi[2*mbs+1] = bdiag[0]+1; /* endof bj and ba array */ 1588 for (i=0; i<mbs; i++) { 1589 /* copy initial fill into linked list */ 1590 nzi = 0; /* nonzeros for active row i */ 1591 nzi = ai[r[i]+1] - ai[r[i]]; 1592 if (!nzi) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i); 1593 nzi_al = adiag[r[i]] - ai[r[i]]; 1594 nzi_au = ai[r[i]+1] - adiag[r[i]] -1; 1595 /* printf("row %d, nzi_al/au %d %d\n",i,nzi_al,nzi_au); */ 1596 1597 /* load in initial unfactored row */ 1598 ajtmp = aj + ai[r[i]]; 1599 ierr = PetscLLAddPerm(nzi,ajtmp,ic,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1600 ierr = PetscMemzero(rtmp,mbs*bs2*sizeof(PetscScalar));CHKERRQ(ierr); 1601 aatmp = a->a + bs2*ai[r[i]]; 1602 for (j=0; j<nzi; j++) { 1603 ierr = PetscMemcpy(rtmp+bs2*ic[ajtmp[j]],aatmp+bs2*j,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1604 } 1605 1606 /* add pivot rows into linked list */ 1607 row = lnk[mbs]; 1608 while (row < i) { 1609 nzi_bl = bi[row+1] - bi[row] + 1; 1610 bjtmp = bj + bdiag[row+1]+1; /* points to 1st column next to the diagonal in U */ 1611 ierr = PetscLLAddSortedLU(bjtmp,row,nlnk,lnk,lnkbt,i,nzi_bl,im);CHKERRQ(ierr); 1612 nzi += nlnk; 1613 row = lnk[row]; 1614 } 1615 1616 /* copy data from lnk into jtmp, then initialize lnk */ 1617 ierr = PetscLLClean(mbs,mbs,nzi,lnk,jtmp,lnkbt);CHKERRQ(ierr); 1618 1619 /* numerical factorization */ 1620 bjtmp = jtmp; 1621 row = *bjtmp++; /* 1st pivot row */ 1622 1623 while (row < i) { 1624 pc = rtmp + bs2*row; 1625 pv = ba + bs2*bdiag[row]; /* inv(diag) of the pivot row */ 1626 PetscKernel_A_gets_A_times_B(bs,pc,pv,multiplier); /* pc= multiplier = pc*inv(diag[row]) */ 1627 ierr = MatBlockAbs_private(1,bs2,pc,vtmp_abs);CHKERRQ(ierr); 1628 if (vtmp_abs[0] > dt) { /* apply tolerance dropping rule */ 1629 pj = bj + bdiag[row+1] + 1; /* point to 1st entry of U(row,:) */ 1630 pv = ba + bs2*(bdiag[row+1] + 1); 1631 nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries in U(row,:), excluding diagonal */ 1632 for (j=0; j<nz; j++) { 1633 PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); 1634 } 1635 /* ierr = PetscLogFlops(bslog*(nz+1.0)-bs);CHKERRQ(ierr); */ 1636 } 1637 row = *bjtmp++; 1638 } 1639 1640 /* copy sparse rtmp into contiguous vtmp; separate L and U part */ 1641 nzi_bl = 0; j = 0; 1642 while (jtmp[j] < i) { /* L-part. Note: jtmp is sorted */ 1643 ierr = PetscMemcpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 1644 nzi_bl++; j++; 1645 } 1646 nzi_bu = nzi - nzi_bl -1; 1647 /* printf("nzi %d, nzi_bl %d, nzi_bu %d\n",nzi,nzi_bl,nzi_bu); */ 1648 1649 while (j < nzi) { /* U-part */ 1650 ierr = PetscMemcpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 1651 /* 1652 printf(" col %d: ",jtmp[j]); 1653 for (j1=0; j1<bs2; j1++) printf(" %g",*(vtmp+bs2*j+j1)); 1654 printf(" \n"); 1655 */ 1656 j++; 1657 } 1658 1659 ierr = MatBlockAbs_private(nzi,bs2,vtmp,vtmp_abs);CHKERRQ(ierr); 1660 /* 1661 printf(" row %d, nzi %d, vtmp_abs\n",i,nzi); 1662 for (j1=0; j1<nzi; j1++) printf(" (%d %g),",jtmp[j1],vtmp_abs[j1]); 1663 printf(" \n"); 1664 */ 1665 bjtmp = bj + bi[i]; 1666 batmp = ba + bs2*bi[i]; 1667 /* apply level dropping rule to L part */ 1668 ncut = nzi_al + dtcount; 1669 if (ncut < nzi_bl) { 1670 ierr = PetscSortSplitReal(ncut,nzi_bl,vtmp_abs,jtmp);CHKERRQ(ierr); 1671 ierr = PetscSortIntWithScalarArray(ncut,jtmp,vtmp);CHKERRQ(ierr); 1672 } else { 1673 ncut = nzi_bl; 1674 } 1675 for (j=0; j<ncut; j++) { 1676 bjtmp[j] = jtmp[j]; 1677 ierr = PetscMemcpy(batmp+bs2*j,rtmp+bs2*bjtmp[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 1678 /* 1679 printf(" col %d: ",bjtmp[j]); 1680 for (j1=0; j1<bs2; j1++) printf(" %g,",*(batmp+bs2*j+j1)); 1681 printf("\n"); 1682 */ 1683 } 1684 bi[i+1] = bi[i] + ncut; 1685 nzi = ncut + 1; 1686 1687 /* apply level dropping rule to U part */ 1688 ncut = nzi_au + dtcount; 1689 if (ncut < nzi_bu) { 1690 ierr = PetscSortSplitReal(ncut,nzi_bu,vtmp_abs+nzi_bl+1,jtmp+nzi_bl+1);CHKERRQ(ierr); 1691 ierr = PetscSortIntWithScalarArray(ncut,jtmp+nzi_bl+1,vtmp+nzi_bl+1);CHKERRQ(ierr); 1692 } else { 1693 ncut = nzi_bu; 1694 } 1695 nzi += ncut; 1696 1697 /* mark bdiagonal */ 1698 bdiag[i+1] = bdiag[i] - (ncut + 1); 1699 bi[2*mbs - i] = bi[2*mbs - i +1] - (ncut + 1); 1700 1701 bjtmp = bj + bdiag[i]; 1702 batmp = ba + bs2*bdiag[i]; 1703 ierr = PetscMemcpy(batmp,rtmp+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1704 *bjtmp = i; 1705 /* 1706 printf(" diag %d: ",*bjtmp); 1707 for (j=0; j<bs2; j++) { 1708 printf(" %g,",batmp[j]); 1709 } 1710 printf("\n"); 1711 */ 1712 bjtmp = bj + bdiag[i+1]+1; 1713 batmp = ba + (bdiag[i+1]+1)*bs2; 1714 1715 for (k=0; k<ncut; k++) { 1716 bjtmp[k] = jtmp[nzi_bl+1+k]; 1717 ierr = PetscMemcpy(batmp+bs2*k,rtmp+bs2*bjtmp[k],bs2*sizeof(MatScalar));CHKERRQ(ierr); 1718 /* 1719 printf(" col %d:",bjtmp[k]); 1720 for (j1=0; j1<bs2; j1++) printf(" %g,",*(batmp+bs2*k+j1)); 1721 printf("\n"); 1722 */ 1723 } 1724 1725 im[i] = nzi; /* used by PetscLLAddSortedLU() */ 1726 1727 /* invert diagonal block for simplier triangular solves - add shift??? */ 1728 batmp = ba + bs2*bdiag[i]; 1729 ierr = PetscKernel_A_gets_inverse_A(bs,batmp,v_pivots,v_work);CHKERRQ(ierr); 1730 } /* for (i=0; i<mbs; i++) */ 1731 ierr = PetscFree3(v_work,multiplier,v_pivots);CHKERRQ(ierr); 1732 1733 /* printf("end of L %d, beginning of U %d\n",bi[mbs],bdiag[mbs]); */ 1734 if (bi[mbs] >= bdiag[mbs]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"end of L array %d cannot >= the beginning of U array %d",bi[mbs],bdiag[mbs]); 1735 1736 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 1737 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 1738 1739 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 1740 1741 ierr = PetscFree2(im,jtmp);CHKERRQ(ierr); 1742 ierr = PetscFree2(rtmp,vtmp);CHKERRQ(ierr); 1743 1744 ierr = PetscLogFlops(bs2*B->cmap->n);CHKERRQ(ierr); 1745 b->maxnz = b->nz = bi[mbs] + bdiag[0] - bdiag[mbs]; 1746 1747 ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr); 1748 ierr = ISIdentity(isicol,&icol_identity);CHKERRQ(ierr); 1749 if (row_identity && icol_identity) { 1750 B->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering; 1751 } else { 1752 B->ops->solve = MatSolve_SeqBAIJ_N; 1753 } 1754 1755 B->ops->solveadd = 0; 1756 B->ops->solvetranspose = 0; 1757 B->ops->solvetransposeadd = 0; 1758 B->ops->matsolve = 0; 1759 B->assembled = PETSC_TRUE; 1760 B->preallocated = PETSC_TRUE; 1761 PetscFunctionReturn(0); 1762 } 1763