File: | mat/impls/baij/seq/baijfact13.c |
Warning: | line 193, column 9 Value stored to 'ierr' is never read |
[?] Use j/k keys for keyboard navigation
1 | |
2 | /* |
3 | Factorization code for BAIJ format. |
4 | */ |
5 | #include <../src/mat/impls/baij/seq/baij.h> |
6 | #include <petsc/private/kernels/blockinvert.h> |
7 | |
8 | /* |
9 | Version for when blocks are 3 by 3 |
10 | */ |
11 | PetscErrorCode MatLUFactorNumeric_SeqBAIJ_3_inplace(Mat C,Mat A,const MatFactorInfo *info) |
12 | { |
13 | Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data; |
14 | IS isrow = b->row,isicol = b->icol; |
15 | PetscErrorCode ierr; |
16 | const PetscInt *r,*ic; |
17 | PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j; |
18 | PetscInt *ajtmpold,*ajtmp,nz,row,*ai=a->i,*aj=a->j; |
19 | PetscInt *diag_offset = b->diag,idx,*pj; |
20 | MatScalar *pv,*v,*rtmp,*pc,*w,*x; |
21 | MatScalar p1,p2,p3,p4,m1,m2,m3,m4,m5,m6,m7,m8,m9,x1,x2,x3,x4; |
22 | MatScalar p5,p6,p7,p8,p9,x5,x6,x7,x8,x9; |
23 | MatScalar *ba = b->a,*aa = a->a; |
24 | PetscReal shift = info->shiftamount; |
25 | PetscBool allowzeropivot,zeropivotdetected; |
26 | |
27 | PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize < 64)) { petscstack->function[petscstack->currentsize ] = __func__; petscstack->file[petscstack->currentsize] = "/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ; petscstack->line[petscstack->currentsize] = 27; petscstack ->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack ->currentsize++; } if (petscstack) { petscstack->hotdepth += (PETSC_FALSE || petscstack->hotdepth); } ; } while (0) ; ; } while (0); |
28 | ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),28,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
29 | ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),29,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
30 | ierr = PetscMalloc1(9*(n+1),&rtmp)PetscMallocA(1,PETSC_FALSE,30,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,(size_t)(9*(n+1))*sizeof(**(&rtmp)),(&rtmp));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),30,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
31 | allowzeropivot = PetscNot(A->erroriffailure)((A->erroriffailure) ? PETSC_FALSE : PETSC_TRUE); |
32 | |
33 | for (i=0; i<n; i++) { |
34 | nz = bi[i+1] - bi[i]; |
35 | ajtmp = bj + bi[i]; |
36 | for (j=0; j<nz; j++) { |
37 | x = rtmp + 9*ajtmp[j]; |
38 | x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = 0.0; |
39 | } |
40 | /* load in initial (unfactored row) */ |
41 | idx = r[i]; |
42 | nz = ai[idx+1] - ai[idx]; |
43 | ajtmpold = aj + ai[idx]; |
44 | v = aa + 9*ai[idx]; |
45 | for (j=0; j<nz; j++) { |
46 | x = rtmp + 9*ic[ajtmpold[j]]; |
47 | x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3]; |
48 | x[4] = v[4]; x[5] = v[5]; x[6] = v[6]; x[7] = v[7]; x[8] = v[8]; |
49 | v += 9; |
50 | } |
51 | row = *ajtmp++; |
52 | while (row < i) { |
53 | pc = rtmp + 9*row; |
54 | p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3]; |
55 | p5 = pc[4]; p6 = pc[5]; p7 = pc[6]; p8 = pc[7]; p9 = pc[8]; |
56 | if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0 || p5 != 0.0 || |
57 | p6 != 0.0 || p7 != 0.0 || p8 != 0.0 || p9 != 0.0) { |
58 | pv = ba + 9*diag_offset[row]; |
59 | pj = bj + diag_offset[row] + 1; |
60 | x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; |
61 | x5 = pv[4]; x6 = pv[5]; x7 = pv[6]; x8 = pv[7]; x9 = pv[8]; |
62 | pc[0] = m1 = p1*x1 + p4*x2 + p7*x3; |
63 | pc[1] = m2 = p2*x1 + p5*x2 + p8*x3; |
64 | pc[2] = m3 = p3*x1 + p6*x2 + p9*x3; |
65 | |
66 | pc[3] = m4 = p1*x4 + p4*x5 + p7*x6; |
67 | pc[4] = m5 = p2*x4 + p5*x5 + p8*x6; |
68 | pc[5] = m6 = p3*x4 + p6*x5 + p9*x6; |
69 | |
70 | pc[6] = m7 = p1*x7 + p4*x8 + p7*x9; |
71 | pc[7] = m8 = p2*x7 + p5*x8 + p8*x9; |
72 | pc[8] = m9 = p3*x7 + p6*x8 + p9*x9; |
73 | nz = bi[row+1] - diag_offset[row] - 1; |
74 | pv += 9; |
75 | for (j=0; j<nz; j++) { |
76 | x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; |
77 | x5 = pv[4]; x6 = pv[5]; x7 = pv[6]; x8 = pv[7]; x9 = pv[8]; |
78 | x = rtmp + 9*pj[j]; |
79 | x[0] -= m1*x1 + m4*x2 + m7*x3; |
80 | x[1] -= m2*x1 + m5*x2 + m8*x3; |
81 | x[2] -= m3*x1 + m6*x2 + m9*x3; |
82 | |
83 | x[3] -= m1*x4 + m4*x5 + m7*x6; |
84 | x[4] -= m2*x4 + m5*x5 + m8*x6; |
85 | x[5] -= m3*x4 + m6*x5 + m9*x6; |
86 | |
87 | x[6] -= m1*x7 + m4*x8 + m7*x9; |
88 | x[7] -= m2*x7 + m5*x8 + m8*x9; |
89 | x[8] -= m3*x7 + m6*x8 + m9*x9; |
90 | pv += 9; |
91 | } |
92 | ierr = PetscLogFlops(54.0*nz+36.0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),92,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
93 | } |
94 | row = *ajtmp++; |
95 | } |
96 | /* finished row so stick it into b->a */ |
97 | pv = ba + 9*bi[i]; |
98 | pj = bj + bi[i]; |
99 | nz = bi[i+1] - bi[i]; |
100 | for (j=0; j<nz; j++) { |
101 | x = rtmp + 9*pj[j]; |
102 | pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3]; |
103 | pv[4] = x[4]; pv[5] = x[5]; pv[6] = x[6]; pv[7] = x[7]; pv[8] = x[8]; |
104 | pv += 9; |
105 | } |
106 | /* invert diagonal block */ |
107 | w = ba + 9*diag_offset[i]; |
108 | ierr = PetscKernel_A_gets_inverse_A_3(w,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),108,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
109 | if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; |
110 | } |
111 | |
112 | ierr = PetscFree(rtmp)((*PetscTrFree)((void*)(rtmp),112,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ) || ((rtmp) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),112,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
113 | ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),113,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
114 | ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),114,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
115 | |
116 | C->ops->solve = MatSolve_SeqBAIJ_3_inplace; |
117 | C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3_inplace; |
118 | C->assembled = PETSC_TRUE; |
119 | |
120 | ierr = PetscLogFlops(1.333333333333*3*3*3*b->mbs);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),120,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */ |
121 | 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); |
122 | } |
123 | |
124 | /* MatLUFactorNumeric_SeqBAIJ_3 - |
125 | copied from MatLUFactorNumeric_SeqBAIJ_N_inplace() and manually re-implemented |
126 | PetscKernel_A_gets_A_times_B() |
127 | PetscKernel_A_gets_A_minus_B_times_C() |
128 | PetscKernel_A_gets_inverse_A() |
129 | */ |
130 | PetscErrorCode MatLUFactorNumeric_SeqBAIJ_3(Mat B,Mat A,const MatFactorInfo *info) |
131 | { |
132 | Mat C =B; |
133 | Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data; |
134 | IS isrow = b->row,isicol = b->icol; |
135 | PetscErrorCode ierr; |
136 | const PetscInt *r,*ic; |
137 | PetscInt i,j,k,nz,nzL,row; |
138 | const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j; |
139 | const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag,*pj,bs2=a->bs2; |
140 | MatScalar *rtmp,*pc,*mwork,*v,*pv,*aa=a->a; |
141 | PetscInt flg; |
142 | PetscReal shift = info->shiftamount; |
143 | PetscBool allowzeropivot,zeropivotdetected; |
144 | |
145 | PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize < 64)) { petscstack->function[petscstack->currentsize ] = __func__; petscstack->file[petscstack->currentsize] = "/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ; petscstack->line[petscstack->currentsize] = 145; petscstack ->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack ->currentsize++; } if (petscstack) { petscstack->hotdepth += (PETSC_FALSE || petscstack->hotdepth); } ; } while (0) ; ; } while (0); |
146 | ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),146,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
147 | ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),147,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
148 | allowzeropivot = PetscNot(A->erroriffailure)((A->erroriffailure) ? PETSC_FALSE : PETSC_TRUE); |
149 | |
150 | /* generate work space needed by the factorization */ |
151 | ierr = PetscMalloc2(bs2*n,&rtmp,bs2,&mwork)PetscMallocA(2,PETSC_FALSE,151,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.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),151,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
152 | ierr = PetscArrayzero(rtmp,bs2*n)PetscMemzero(rtmp,(bs2*n)*sizeof(*(rtmp)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),152,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
153 | |
154 | for (i=0; i<n; i++) { |
155 | /* zero rtmp */ |
156 | /* L part */ |
157 | nz = bi[i+1] - bi[i]; |
158 | bjtmp = bj + bi[i]; |
159 | for (j=0; j<nz; j++) { |
160 | 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),160,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
161 | } |
162 | |
163 | /* U part */ |
164 | nz = bdiag[i] - bdiag[i+1]; |
165 | bjtmp = bj + bdiag[i+1]+1; |
166 | for (j=0; j<nz; j++) { |
167 | 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),167,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
168 | } |
169 | |
170 | /* load in initial (unfactored row) */ |
171 | nz = ai[r[i]+1] - ai[r[i]]; |
172 | ajtmp = aj + ai[r[i]]; |
173 | v = aa + bs2*ai[r[i]]; |
174 | for (j=0; j<nz; j++) { |
175 | 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),175,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
176 | } |
177 | |
178 | /* elimination */ |
179 | bjtmp = bj + bi[i]; |
180 | nzL = bi[i+1] - bi[i]; |
181 | for (k = 0; k < nzL; k++) { |
182 | row = bjtmp[k]; |
183 | pc = rtmp + bs2*row; |
184 | for (flg=0,j=0; j<bs2; j++) { |
185 | if (pc[j]!=0.0) { |
186 | flg = 1; |
187 | break; |
188 | } |
189 | } |
190 | if (flg) { |
191 | pv = b->a + bs2*bdiag[row]; |
192 | /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ |
193 | ierr = PetscKernel_A_gets_A_times_B_3(pc,pv,mwork)0; { ierr = ((sizeof(*(mwork)) != sizeof(*(pc))) || PetscMemcpy (mwork,pc,(9)*sizeof(*(mwork))));do {if (__builtin_expect(!!( ierr),0)) return PetscError(((MPI_Comm)0x44000001),193,__func__ ,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); pc[0] = mwork[0]*pv [0] + mwork[3]*pv[1] + mwork[6]*pv[2]; pc[1] = mwork[1]*pv[0] + mwork[4]*pv[1] + mwork[7]*pv[2]; pc[2] = mwork[2]*pv[0] + mwork [5]*pv[1] + mwork[8]*pv[2]; pc[3] = mwork[0]*pv[3] + mwork[3] *pv[4] + mwork[6]*pv[5]; pc[4] = mwork[1]*pv[3] + mwork[4]*pv [4] + mwork[7]*pv[5]; pc[5] = mwork[2]*pv[3] + mwork[5]*pv[4] + mwork[8]*pv[5]; pc[6] = mwork[0]*pv[6] + mwork[3]*pv[7] + mwork [6]*pv[8]; pc[7] = mwork[1]*pv[6] + mwork[4]*pv[7] + mwork[7] *pv[8]; pc[8] = mwork[2]*pv[6] + mwork[5]*pv[7] + mwork[8]*pv [8]; };CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),193,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
Value stored to 'ierr' is never read | |
194 | |
195 | pj = b->j + bdiag[row+1] + 1; /* beginning of U(row,:) */ |
196 | pv = b->a + bs2*(bdiag[row+1]+1); |
197 | nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries in U(row,:) excluding diag */ |
198 | for (j=0; j<nz; j++) { |
199 | /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ |
200 | /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ |
201 | v = rtmp + bs2*pj[j]; |
202 | ierr = PetscKernel_A_gets_A_minus_B_times_C_3(v,pc,pv)0; { v[0] -= pc[0]*pv[0] + pc[3]*pv[1] + pc[6]*pv[2]; v[1] -= pc[1]*pv[0] + pc[4]*pv[1] + pc[7]*pv[2]; v[2] -= pc[2]*pv[0] + pc[5]*pv[1] + pc[8]*pv[2]; v[3] -= pc[0]*pv[3] + pc[3]*pv[ 4] + pc[6]*pv[5]; v[4] -= pc[1]*pv[3] + pc[4]*pv[4] + pc[7]*pv [5]; v[5] -= pc[2]*pv[3] + pc[5]*pv[4] + pc[8]*pv[5]; v[6] -= pc[0]*pv[6] + pc[3]*pv[7] + pc[6]*pv[8]; v[7] -= pc[1]*pv[6] + pc[4]*pv[7] + pc[7]*pv[8]; v[8] -= pc[2]*pv[6] + pc[5]*pv[ 7] + pc[8]*pv[8]; };CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),202,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
203 | pv += bs2; |
204 | } |
205 | ierr = PetscLogFlops(54*nz+45);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),205,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ |
206 | } |
207 | } |
208 | |
209 | /* finished row so stick it into b->a */ |
210 | /* L part */ |
211 | pv = b->a + bs2*bi[i]; |
212 | pj = b->j + bi[i]; |
213 | nz = bi[i+1] - bi[i]; |
214 | for (j=0; j<nz; j++) { |
215 | 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),215,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
216 | } |
217 | |
218 | /* Mark diagonal and invert diagonal for simplier triangular solves */ |
219 | pv = b->a + bs2*bdiag[i]; |
220 | pj = b->j + bdiag[i]; |
221 | 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),221,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
222 | ierr = PetscKernel_A_gets_inverse_A_3(pv,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),222,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
223 | if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; |
224 | |
225 | /* U part */ |
226 | pj = b->j + bdiag[i+1] + 1; |
227 | pv = b->a + bs2*(bdiag[i+1]+1); |
228 | nz = bdiag[i] - bdiag[i+1] - 1; |
229 | for (j=0; j<nz; j++) { |
230 | 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),230,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
231 | } |
232 | } |
233 | |
234 | ierr = PetscFree2(rtmp,mwork)PetscFreeA(2,234,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,&(rtmp),&(mwork));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),234,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
235 | ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),235,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
236 | ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),236,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
237 | |
238 | C->ops->solve = MatSolve_SeqBAIJ_3; |
239 | C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3; |
240 | C->assembled = PETSC_TRUE; |
241 | |
242 | ierr = PetscLogFlops(1.333333333333*3*3*3*n);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),242,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */ |
243 | 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); |
244 | } |
245 | |
246 | PetscErrorCode MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info) |
247 | { |
248 | Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data; |
249 | PetscErrorCode ierr; |
250 | PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j; |
251 | PetscInt *ajtmpold,*ajtmp,nz,row; |
252 | PetscInt *diag_offset = b->diag,*ai=a->i,*aj=a->j,*pj; |
253 | MatScalar *pv,*v,*rtmp,*pc,*w,*x; |
254 | MatScalar p1,p2,p3,p4,m1,m2,m3,m4,m5,m6,m7,m8,m9,x1,x2,x3,x4; |
255 | MatScalar p5,p6,p7,p8,p9,x5,x6,x7,x8,x9; |
256 | MatScalar *ba = b->a,*aa = a->a; |
257 | PetscReal shift = info->shiftamount; |
258 | PetscBool allowzeropivot,zeropivotdetected; |
259 | |
260 | PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize < 64)) { petscstack->function[petscstack->currentsize ] = __func__; petscstack->file[petscstack->currentsize] = "/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ; petscstack->line[petscstack->currentsize] = 260; petscstack ->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack ->currentsize++; } if (petscstack) { petscstack->hotdepth += (PETSC_FALSE || petscstack->hotdepth); } ; } while (0) ; ; } while (0); |
261 | ierr = PetscMalloc1(9*(n+1),&rtmp)PetscMallocA(1,PETSC_FALSE,261,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,(size_t)(9*(n+1))*sizeof(**(&rtmp)),(&rtmp));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),261,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
262 | allowzeropivot = PetscNot(A->erroriffailure)((A->erroriffailure) ? PETSC_FALSE : PETSC_TRUE); |
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+9*ajtmp[j]; |
269 | x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = 0.0; |
270 | } |
271 | /* load in initial (unfactored row) */ |
272 | nz = ai[i+1] - ai[i]; |
273 | ajtmpold = aj + ai[i]; |
274 | v = aa + 9*ai[i]; |
275 | for (j=0; j<nz; j++) { |
276 | x = rtmp+9*ajtmpold[j]; |
277 | x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3]; |
278 | x[4] = v[4]; x[5] = v[5]; x[6] = v[6]; x[7] = v[7]; x[8] = v[8]; |
279 | v += 9; |
280 | } |
281 | row = *ajtmp++; |
282 | while (row < i) { |
283 | pc = rtmp + 9*row; |
284 | p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3]; |
285 | p5 = pc[4]; p6 = pc[5]; p7 = pc[6]; p8 = pc[7]; p9 = pc[8]; |
286 | if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0 || p5 != 0.0 || |
287 | p6 != 0.0 || p7 != 0.0 || p8 != 0.0 || p9 != 0.0) { |
288 | pv = ba + 9*diag_offset[row]; |
289 | pj = bj + diag_offset[row] + 1; |
290 | x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; |
291 | x5 = pv[4]; x6 = pv[5]; x7 = pv[6]; x8 = pv[7]; x9 = pv[8]; |
292 | pc[0] = m1 = p1*x1 + p4*x2 + p7*x3; |
293 | pc[1] = m2 = p2*x1 + p5*x2 + p8*x3; |
294 | pc[2] = m3 = p3*x1 + p6*x2 + p9*x3; |
295 | |
296 | pc[3] = m4 = p1*x4 + p4*x5 + p7*x6; |
297 | pc[4] = m5 = p2*x4 + p5*x5 + p8*x6; |
298 | pc[5] = m6 = p3*x4 + p6*x5 + p9*x6; |
299 | |
300 | pc[6] = m7 = p1*x7 + p4*x8 + p7*x9; |
301 | pc[7] = m8 = p2*x7 + p5*x8 + p8*x9; |
302 | pc[8] = m9 = p3*x7 + p6*x8 + p9*x9; |
303 | |
304 | nz = bi[row+1] - diag_offset[row] - 1; |
305 | pv += 9; |
306 | for (j=0; j<nz; j++) { |
307 | x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; |
308 | x5 = pv[4]; x6 = pv[5]; x7 = pv[6]; x8 = pv[7]; x9 = pv[8]; |
309 | x = rtmp + 9*pj[j]; |
310 | x[0] -= m1*x1 + m4*x2 + m7*x3; |
311 | x[1] -= m2*x1 + m5*x2 + m8*x3; |
312 | x[2] -= m3*x1 + m6*x2 + m9*x3; |
313 | |
314 | x[3] -= m1*x4 + m4*x5 + m7*x6; |
315 | x[4] -= m2*x4 + m5*x5 + m8*x6; |
316 | x[5] -= m3*x4 + m6*x5 + m9*x6; |
317 | |
318 | x[6] -= m1*x7 + m4*x8 + m7*x9; |
319 | x[7] -= m2*x7 + m5*x8 + m8*x9; |
320 | x[8] -= m3*x7 + m6*x8 + m9*x9; |
321 | pv += 9; |
322 | } |
323 | ierr = PetscLogFlops(54.0*nz+36.0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),323,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
324 | } |
325 | row = *ajtmp++; |
326 | } |
327 | /* finished row so stick it into b->a */ |
328 | pv = ba + 9*bi[i]; |
329 | pj = bj + bi[i]; |
330 | nz = bi[i+1] - bi[i]; |
331 | for (j=0; j<nz; j++) { |
332 | x = rtmp+9*pj[j]; |
333 | pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3]; |
334 | pv[4] = x[4]; pv[5] = x[5]; pv[6] = x[6]; pv[7] = x[7]; pv[8] = x[8]; |
335 | pv += 9; |
336 | } |
337 | /* invert diagonal block */ |
338 | w = ba + 9*diag_offset[i]; |
339 | ierr = PetscKernel_A_gets_inverse_A_3(w,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),339,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
340 | if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; |
341 | } |
342 | |
343 | ierr = PetscFree(rtmp)((*PetscTrFree)((void*)(rtmp),343,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ) || ((rtmp) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),343,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
344 | |
345 | C->ops->solve = MatSolve_SeqBAIJ_3_NaturalOrdering_inplace; |
346 | C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3_NaturalOrdering_inplace; |
347 | C->assembled = PETSC_TRUE; |
348 | |
349 | ierr = PetscLogFlops(1.333333333333*3*3*3*b->mbs);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),349,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */ |
350 | 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); |
351 | } |
352 | |
353 | /* |
354 | MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering - |
355 | copied from MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace() |
356 | */ |
357 | PetscErrorCode MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info) |
358 | { |
359 | Mat C =B; |
360 | Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data; |
361 | PetscErrorCode ierr; |
362 | PetscInt i,j,k,nz,nzL,row; |
363 | const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j; |
364 | const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag,*pj,bs2=a->bs2; |
365 | MatScalar *rtmp,*pc,*mwork,*v,*pv,*aa=a->a; |
366 | PetscInt flg; |
367 | PetscReal shift = info->shiftamount; |
368 | PetscBool allowzeropivot,zeropivotdetected; |
369 | |
370 | PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize < 64)) { petscstack->function[petscstack->currentsize ] = __func__; petscstack->file[petscstack->currentsize] = "/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ; petscstack->line[petscstack->currentsize] = 370; petscstack ->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack ->currentsize++; } if (petscstack) { petscstack->hotdepth += (PETSC_FALSE || petscstack->hotdepth); } ; } while (0) ; ; } while (0); |
371 | allowzeropivot = PetscNot(A->erroriffailure)((A->erroriffailure) ? PETSC_FALSE : PETSC_TRUE); |
372 | |
373 | /* generate work space needed by the factorization */ |
374 | ierr = PetscMalloc2(bs2*n,&rtmp,bs2,&mwork)PetscMallocA(2,PETSC_FALSE,374,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.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),374,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
375 | ierr = PetscArrayzero(rtmp,bs2*n)PetscMemzero(rtmp,(bs2*n)*sizeof(*(rtmp)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),375,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
376 | |
377 | for (i=0; i<n; i++) { |
378 | /* zero rtmp */ |
379 | /* L part */ |
380 | nz = bi[i+1] - bi[i]; |
381 | bjtmp = bj + bi[i]; |
382 | for (j=0; j<nz; j++) { |
383 | 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),383,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
384 | } |
385 | |
386 | /* U part */ |
387 | nz = bdiag[i] - bdiag[i+1]; |
388 | bjtmp = bj + bdiag[i+1] + 1; |
389 | for (j=0; j<nz; j++) { |
390 | 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),390,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
391 | } |
392 | |
393 | /* load in initial (unfactored row) */ |
394 | nz = ai[i+1] - ai[i]; |
395 | ajtmp = aj + ai[i]; |
396 | v = aa + bs2*ai[i]; |
397 | for (j=0; j<nz; j++) { |
398 | 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),398,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
399 | } |
400 | |
401 | /* elimination */ |
402 | bjtmp = bj + bi[i]; |
403 | nzL = bi[i+1] - bi[i]; |
404 | for (k=0; k<nzL; k++) { |
405 | row = bjtmp[k]; |
406 | pc = rtmp + bs2*row; |
407 | for (flg=0,j=0; j<bs2; j++) { |
408 | if (pc[j]!=0.0) { |
409 | flg = 1; |
410 | break; |
411 | } |
412 | } |
413 | if (flg) { |
414 | pv = b->a + bs2*bdiag[row]; |
415 | /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ |
416 | ierr = PetscKernel_A_gets_A_times_B_3(pc,pv,mwork)0; { ierr = ((sizeof(*(mwork)) != sizeof(*(pc))) || PetscMemcpy (mwork,pc,(9)*sizeof(*(mwork))));do {if (__builtin_expect(!!( ierr),0)) return PetscError(((MPI_Comm)0x44000001),416,__func__ ,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); pc[0] = mwork[0]*pv [0] + mwork[3]*pv[1] + mwork[6]*pv[2]; pc[1] = mwork[1]*pv[0] + mwork[4]*pv[1] + mwork[7]*pv[2]; pc[2] = mwork[2]*pv[0] + mwork [5]*pv[1] + mwork[8]*pv[2]; pc[3] = mwork[0]*pv[3] + mwork[3] *pv[4] + mwork[6]*pv[5]; pc[4] = mwork[1]*pv[3] + mwork[4]*pv [4] + mwork[7]*pv[5]; pc[5] = mwork[2]*pv[3] + mwork[5]*pv[4] + mwork[8]*pv[5]; pc[6] = mwork[0]*pv[6] + mwork[3]*pv[7] + mwork [6]*pv[8]; pc[7] = mwork[1]*pv[6] + mwork[4]*pv[7] + mwork[7] *pv[8]; pc[8] = mwork[2]*pv[6] + mwork[5]*pv[7] + mwork[8]*pv [8]; };CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),416,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
417 | |
418 | pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */ |
419 | pv = b->a + bs2*(bdiag[row+1]+1); |
420 | nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries in U(row,:) excluding diag */ |
421 | for (j=0; j<nz; j++) { |
422 | /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ |
423 | /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ |
424 | v = rtmp + bs2*pj[j]; |
425 | ierr = PetscKernel_A_gets_A_minus_B_times_C_3(v,pc,pv)0; { v[0] -= pc[0]*pv[0] + pc[3]*pv[1] + pc[6]*pv[2]; v[1] -= pc[1]*pv[0] + pc[4]*pv[1] + pc[7]*pv[2]; v[2] -= pc[2]*pv[0] + pc[5]*pv[1] + pc[8]*pv[2]; v[3] -= pc[0]*pv[3] + pc[3]*pv[ 4] + pc[6]*pv[5]; v[4] -= pc[1]*pv[3] + pc[4]*pv[4] + pc[7]*pv [5]; v[5] -= pc[2]*pv[3] + pc[5]*pv[4] + pc[8]*pv[5]; v[6] -= pc[0]*pv[6] + pc[3]*pv[7] + pc[6]*pv[8]; v[7] -= pc[1]*pv[6] + pc[4]*pv[7] + pc[7]*pv[8]; v[8] -= pc[2]*pv[6] + pc[5]*pv[ 7] + pc[8]*pv[8]; };CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),425,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
426 | pv += bs2; |
427 | } |
428 | ierr = PetscLogFlops(54*nz+45);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),428,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ |
429 | } |
430 | } |
431 | |
432 | /* finished row so stick it into b->a */ |
433 | /* L part */ |
434 | pv = b->a + bs2*bi[i]; |
435 | pj = b->j + bi[i]; |
436 | nz = bi[i+1] - bi[i]; |
437 | for (j=0; j<nz; j++) { |
438 | 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),438,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
439 | } |
440 | |
441 | /* Mark diagonal and invert diagonal for simplier triangular solves */ |
442 | pv = b->a + bs2*bdiag[i]; |
443 | pj = b->j + bdiag[i]; |
444 | 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),444,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
445 | ierr = PetscKernel_A_gets_inverse_A_3(pv,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),445,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
446 | if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; |
447 | |
448 | /* U part */ |
449 | pv = b->a + bs2*(bdiag[i+1]+1); |
450 | pj = b->j + bdiag[i+1]+1; |
451 | nz = bdiag[i] - bdiag[i+1] - 1; |
452 | for (j=0; j<nz; j++) { |
453 | 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),453,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
454 | } |
455 | } |
456 | ierr = PetscFree2(rtmp,mwork)PetscFreeA(2,456,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,&(rtmp),&(mwork));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),456,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
457 | |
458 | C->ops->solve = MatSolve_SeqBAIJ_3_NaturalOrdering; |
459 | C->ops->forwardsolve = MatForwardSolve_SeqBAIJ_3_NaturalOrdering; |
460 | C->ops->backwardsolve = MatBackwardSolve_SeqBAIJ_3_NaturalOrdering; |
461 | C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_3_NaturalOrdering; |
462 | C->assembled = PETSC_TRUE; |
463 | |
464 | ierr = PetscLogFlops(1.333333333333*3*3*3*n);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),464,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/baij/seq/baijfact13.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */ |
465 | 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); |
466 | } |
467 |