Bug Summary

File:mat/impls/baij/seq/baijfact.c
Warning:line 85, column 28
The result of the left shift is undefined because the right operand is negative

Annotated Source Code

[?] Use j/k keys for keyboard navigation

/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact.c

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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),72,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
73
74 pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */
75 pv = b->a + bs2*(bdiag[row+1]+1);
76 nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */
77 for (j=0; j<nz; j++) {
78 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
79 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
80 v = rtmp + 4*pj[j];
81 ierr = PetscKernel_A_gets_A_minus_B_times_C_2(v,pc,pv);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),81,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),185,__func__,"/sandbox/petsc/petsc.master/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);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),194,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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)fabs(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])fabs(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.master/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.master/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.master/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.master/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.master/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])fabs(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])fabs(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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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])fabs(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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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])fabs(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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
,73,PETSC_ERROR_INITIAL,"Matrix is missing diagonal entry %D"
,i)
;
1
Assuming 'missing' is 0
2
Taking false branch
984
985 if (bs > 1) {
3
Assuming 'bs' is <= 1
4
Taking false branch
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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1033 for (i=0; i<am; i++) {
5
Assuming 'i' is < 'am'
6
Loop condition is true. Entering loop body
7
Assuming 'i' is >= 'am'
8
Loop condition is false. Execution continues on line 1038
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.master/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.master/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.master/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.master/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 */
9
Loop condition is true. Entering loop body
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++) {
10
Assuming 'j' is < 'ncols'
11
Loop condition is true. Entering loop body
14
Assuming 'j' is < 'ncols'
15
Loop condition is true. Entering loop body
18
Assuming 'j' is >= 'ncols'
19
Loop condition is false. Execution continues on line 1063
1056 i = rip[*(aj + ai[rip[k]] + j)];
1057 if (i >= k) { /* only take upper triangular entry */
12
Assuming 'i' is >= 'k'
13
Taking true branch
16
Assuming 'i' is >= 'k'
17
Taking true branch
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.master/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
20
Within the expansion of the macro 'PetscIncompleteLLAdd':
a
Calling 'PetscBTLookupSet'
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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),1393,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),1406,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),1409,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),1448,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),1461,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),1464,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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])fabs(blockarray[i*nbs+j])) absarray[i] = PetscAbsScalar(blockarray[i*nbs+j])fabs(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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),1628,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/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.master/src/mat/impls/baij/seq/baijfact.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),1635,__func__,"/sandbox/petsc/petsc.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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.master/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}

/sandbox/petsc/petsc.master/include/petscbt.h

1
2#if !defined(PETSCBT_H)
3#define PETSCBT_H
4
5#include <petscconf.h>
6#include <petscviewer.h>
7
8/*S
9 PetscBT - PETSc bitarrays
10
11 Level: advanced
12
13 PetscBTCreate(m,&bt) - creates a bit array with enough room to hold m values
14 PetscBTDestroy(&bt) - destroys the bit array
15 PetscBTMemzero(m,bt) - zeros the entire bit array (sets all values to false)
16 PetscBTSet(bt,index) - sets a particular entry as true
17 PetscBTClear(bt,index) - sets a particular entry as false
18 PetscBTLookup(bt,index) - returns the value
19 PetscBTLookupSet(bt,index) - returns the value and then sets it true
20 PetscBTLookupClear(bt,index) - returns the value and then sets it false
21 PetscBTLength(m) - returns number of bytes in array with m bits
22 PetscBTView(m,bt,viewer) - prints all the entries in a bit array
23
24 We do not currently check error flags on PetscBTSet(), PetscBTClear(), PetscBTLookup(),
25 PetcBTLookupSet(), PetscBTLength() cause error checking would cost hundreds more cycles then
26 the operation.
27
28S*/
29typedef char* PetscBT;
30
31
32PETSC_STATIC_INLINEstatic inline PetscInt PetscBTLength(PetscInt m)
33{
34 return m/PETSC_BITS_PER_BYTE8+1;
35}
36
37PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTMemzero(PetscInt m,PetscBT array)
38{
39 return PetscMemzero(array,sizeof(char)*((size_t)m/PETSC_BITS_PER_BYTE8+1));
40}
41
42PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTDestroy(PetscBT *array)
43{
44 return PetscFree(*array)((*PetscTrFree)((void*)(*array),44,__func__,"/sandbox/petsc/petsc.master/include/petscbt.h"
) || ((*array) = 0,0))
;
45}
46
47PETSC_STATIC_INLINEstatic inline char PetscBTLookup(PetscBT array,PetscInt index)
48{
49 char BT_mask,BT_c;
50 PetscInt BT_idx;
51
52 BT_idx = index/PETSC_BITS_PER_BYTE8;
53 BT_c = array[BT_idx];
54 BT_mask = (char)(1 << index%PETSC_BITS_PER_BYTE8);
55 return (char)(BT_c & BT_mask);
56}
57
58PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTView(PetscInt m,const PetscBT bt,PetscViewer viewer)
59{
60 PetscInt i;
61 PetscErrorCode ierr;
62
63 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PETSC_COMM_SELF((MPI_Comm)0x44000001),&viewer);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),63,__func__,"/sandbox/petsc/petsc.master/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;}
64 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),64,__func__,"/sandbox/petsc/petsc.master/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
65 for (i=0; i<m; i++) {
66 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%D %d\n",i,(int)PetscBTLookup(bt,i));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),66,__func__,"/sandbox/petsc/petsc.master/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
67 }
68 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),68,__func__,"/sandbox/petsc/petsc.master/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
69 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),69,__func__,"/sandbox/petsc/petsc.master/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
70 return 0;
71}
72
73PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTCreate(PetscInt m,PetscBT *array)
74{
75 return PetscMalloc1((size_t)m/PETSC_BITS_PER_BYTE+1,array)PetscMallocA(1,PETSC_FALSE,75,__func__,"/sandbox/petsc/petsc.master/include/petscbt.h"
,(size_t)((size_t)m/8 +1)*sizeof(**(array)),(array))
|| PetscBTMemzero(m,*array);
76}
77
78PETSC_STATIC_INLINEstatic inline char PetscBTLookupSet(PetscBT array,PetscInt index)
79{
80 char BT_mask,BT_c;
81 PetscInt BT_idx;
82
83 BT_idx = index/PETSC_BITS_PER_BYTE8;
84 BT_c = array[BT_idx];
85 BT_mask = (char)(1 << index%PETSC_BITS_PER_BYTE8);
21
The result of the left shift is undefined because the right operand is negative
86 array[BT_idx] = (char)(BT_c | BT_mask);
87 return (char)(BT_c & BT_mask);
88}
89
90PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTSet(PetscBT array,PetscInt index)
91{
92 char BT_mask,BT_c;
93 PetscInt BT_idx;
94
95 BT_idx = index/PETSC_BITS_PER_BYTE8;
96 BT_c = array[BT_idx];
97 BT_mask = (char)(1 << index%PETSC_BITS_PER_BYTE8);
98 array[BT_idx] = (char)(BT_c | BT_mask);
99 return 0;
100}
101
102PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTNegate(PetscBT array,PetscInt index)
103{
104 const PetscInt BT_idx = index/PETSC_BITS_PER_BYTE8;
105 const char BT_mask = (char)(1 << index%PETSC_BITS_PER_BYTE8);
106
107 array[BT_idx] ^= BT_mask;
108 return 0;
109}
110
111PETSC_STATIC_INLINEstatic inline char PetscBTLookupClear(PetscBT array,PetscInt index)
112{
113 char BT_mask,BT_c;
114 PetscInt BT_idx;
115
116 BT_idx = index/PETSC_BITS_PER_BYTE8;
117 BT_c = array[BT_idx];
118 BT_mask = (char)(1 << index%PETSC_BITS_PER_BYTE8);
119 array[BT_idx] = (char)(BT_c & ~BT_mask);
120 return (char)(BT_c & BT_mask);
121}
122
123PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTClear(PetscBT array,PetscInt index)
124{
125 char BT_mask,BT_c;
126 PetscInt BT_idx;
127
128 BT_idx = index/PETSC_BITS_PER_BYTE8;
129 BT_c = array[BT_idx];
130 BT_mask = (char)(1 << index%PETSC_BITS_PER_BYTE8);
131 array[BT_idx] = (char)(BT_c & ~BT_mask);
132 return 0;
133}
134
135
136#endif