Bug Summary

File:mat/impls/baij/seq/baijfact.c
Warning:line 72, column 9
Value stored to 'ierr' is never read

Annotated Source Code

[?] Use j/k keys for keyboard navigation

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