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