Actual source code: baijfact4.c
2: /*
3: Factorization code for BAIJ format.
4: */
5: #include <../src/mat/impls/baij/seq/baij.h>
6: #include <petsc/private/kernels/blockinvert.h>
8: /* ----------------------------------------------------------- */
9: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N_inplace(Mat C,Mat A,const MatFactorInfo *info)
10: {
11: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data;
12: IS isrow = b->row,isicol = b->icol;
13: const PetscInt *r,*ic;
14: PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j;
15: PetscInt *ajtmpold,*ajtmp,nz,row,*ai=a->i,*aj=a->j,k,flg;
16: PetscInt *diag_offset=b->diag,diag,bs=A->rmap->bs,bs2 = a->bs2,*pj,*v_pivots;
17: MatScalar *ba = b->a,*aa = a->a,*pv,*v,*rtmp,*multiplier,*v_work,*pc,*w;
18: PetscBool allowzeropivot,zeropivotdetected;
20: ISGetIndices(isrow,&r);
21: ISGetIndices(isicol,&ic);
22: allowzeropivot = PetscNot(A->erroriffailure);
24: PetscCalloc1(bs2*(n+1),&rtmp);
25: /* generate work space needed by dense LU factorization */
26: PetscMalloc3(bs,&v_work,bs2,&multiplier,bs,&v_pivots);
28: for (i=0; i<n; i++) {
29: nz = bi[i+1] - bi[i];
30: ajtmp = bj + bi[i];
31: for (j=0; j<nz; j++) {
32: PetscArrayzero(rtmp+bs2*ajtmp[j],bs2);
33: }
34: /* load in initial (unfactored row) */
35: nz = ai[r[i]+1] - ai[r[i]];
36: ajtmpold = aj + ai[r[i]];
37: v = aa + bs2*ai[r[i]];
38: for (j=0; j<nz; j++) {
39: PetscArraycpy(rtmp+bs2*ic[ajtmpold[j]],v+bs2*j,bs2);
40: }
41: row = *ajtmp++;
42: while (row < i) {
43: pc = rtmp + bs2*row;
44: /* if (*pc) { */
45: for (flg=0,k=0; k<bs2; k++) {
46: if (pc[k]!=0.0) {
47: flg = 1;
48: break;
49: }
50: }
51: if (flg) {
52: pv = ba + bs2*diag_offset[row];
53: pj = bj + diag_offset[row] + 1;
54: PetscKernel_A_gets_A_times_B(bs,pc,pv,multiplier);
55: nz = bi[row+1] - diag_offset[row] - 1;
56: pv += bs2;
57: for (j=0; j<nz; j++) {
58: PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j);
59: }
60: PetscLogFlops(2.0*bs*bs2*(nz+1.0)-bs);
61: }
62: row = *ajtmp++;
63: }
64: /* finished row so stick it into b->a */
65: pv = ba + bs2*bi[i];
66: pj = bj + bi[i];
67: nz = bi[i+1] - bi[i];
68: for (j=0; j<nz; j++) {
69: PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);
70: }
71: diag = diag_offset[i] - bi[i];
72: /* invert diagonal block */
73: w = pv + bs2*diag;
75: PetscKernel_A_gets_inverse_A(bs,w,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
76: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
77: }
79: PetscFree(rtmp);
80: PetscFree3(v_work,multiplier,v_pivots);
81: ISRestoreIndices(isicol,&ic);
82: ISRestoreIndices(isrow,&r);
84: C->ops->solve = MatSolve_SeqBAIJ_N_inplace;
85: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N_inplace;
86: C->assembled = PETSC_TRUE;
88: PetscLogFlops(1.333333333333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
89: return 0;
90: }