File: | mat/impls/baij/seq/baijfact9.c |
Warning: | line 277, 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 | /* |
10 | Version for when blocks are 5 by 5 |
11 | */ |
12 | PetscErrorCode MatLUFactorNumeric_SeqBAIJ_5_inplace(Mat C,Mat A,const MatFactorInfo *info) |
13 | { |
14 | Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data; |
15 | IS isrow = b->row,isicol = b->icol; |
16 | PetscErrorCode ierr; |
17 | const PetscInt *r,*ic,*bi = b->i,*bj = b->j,*ajtmpold,*ajtmp; |
18 | PetscInt i,j,n = a->mbs,nz,row,idx,ipvt[5]; |
19 | const PetscInt *diag_offset = b->diag,*ai=a->i,*aj=a->j,*pj; |
20 | MatScalar *w,*pv,*rtmp,*x,*pc; |
21 | const MatScalar *v,*aa = a->a; |
22 | MatScalar p1,p2,p3,p4,m1,m2,m3,m4,m5,m6,m7,m8,m9,x1,x2,x3,x4; |
23 | MatScalar p5,p6,p7,p8,p9,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16; |
24 | MatScalar x17,x18,x19,x20,x21,x22,x23,x24,x25,p10,p11,p12,p13,p14; |
25 | MatScalar p15,p16,p17,p18,p19,p20,p21,p22,p23,p24,p25,m10,m11,m12; |
26 | MatScalar m13,m14,m15,m16,m17,m18,m19,m20,m21,m22,m23,m24,m25; |
27 | MatScalar *ba = b->a,work[25]; |
28 | PetscReal shift = info->shiftamount; |
29 | PetscBool allowzeropivot,zeropivotdetected; |
30 | |
31 | 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/baijfact9.c" ; petscstack->line[petscstack->currentsize] = 31; petscstack ->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack ->currentsize++; } if (petscstack) { petscstack->hotdepth += (PETSC_FALSE || petscstack->hotdepth); } ; } while (0) ; ; } while (0); |
32 | allowzeropivot = PetscNot(A->erroriffailure)((A->erroriffailure) ? PETSC_FALSE : PETSC_TRUE); |
33 | ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),33,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
34 | ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),34,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
35 | ierr = PetscMalloc1(25*(n+1),&rtmp)PetscMallocA(1,PETSC_FALSE,35,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,(size_t)(25*(n+1))*sizeof(**(&rtmp)),(&rtmp));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),35,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
36 | |
37 | #define PETSC_USE_MEMZERO1 1 |
38 | #define PETSC_USE_MEMCPY1 1 |
39 | |
40 | for (i=0; i<n; i++) { |
41 | nz = bi[i+1] - bi[i]; |
42 | ajtmp = bj + bi[i]; |
43 | for (j=0; j<nz; j++) { |
44 | #if defined(PETSC_USE_MEMZERO1) |
45 | ierr = PetscArrayzero(rtmp+25*ajtmp[j],25)PetscMemzero(rtmp+25*ajtmp[j],(25)*sizeof(*(rtmp+25*ajtmp[j]) ));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),45,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
46 | #else |
47 | x = rtmp+25*ajtmp[j]; |
48 | x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = x[9] = 0.0; |
49 | x[10] = x[11] = x[12] = x[13] = x[14] = x[15] = x[16] = x[17] = 0.0; |
50 | x[18] = x[19] = x[20] = x[21] = x[22] = x[23] = x[24] = 0.0; |
51 | #endif |
52 | } |
53 | /* load in initial (unfactored row) */ |
54 | idx = r[i]; |
55 | nz = ai[idx+1] - ai[idx]; |
56 | ajtmpold = aj + ai[idx]; |
57 | v = aa + 25*ai[idx]; |
58 | for (j=0; j<nz; j++) { |
59 | #if defined(PETSC_USE_MEMCPY1) |
60 | ierr = PetscArraycpy(rtmp+25*ic[ajtmpold[j]],v,25)((sizeof(*(rtmp+25*ic[ajtmpold[j]])) != sizeof(*(v))) || PetscMemcpy (rtmp+25*ic[ajtmpold[j]],v,(25)*sizeof(*(rtmp+25*ic[ajtmpold[ j]]))));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),60,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
61 | #else |
62 | x = rtmp+25*ic[ajtmpold[j]]; |
63 | x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3]; |
64 | x[4] = v[4]; x[5] = v[5]; x[6] = v[6]; x[7] = v[7]; x[8] = v[8]; |
65 | x[9] = v[9]; x[10] = v[10]; x[11] = v[11]; x[12] = v[12]; x[13] = v[13]; |
66 | x[14] = v[14]; x[15] = v[15]; x[16] = v[16]; x[17] = v[17]; |
67 | x[18] = v[18]; x[19] = v[19]; x[20] = v[20]; x[21] = v[21]; |
68 | x[22] = v[22]; x[23] = v[23]; x[24] = v[24]; |
69 | #endif |
70 | v += 25; |
71 | } |
72 | row = *ajtmp++; |
73 | while (row < i) { |
74 | pc = rtmp + 25*row; |
75 | p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3]; |
76 | p5 = pc[4]; p6 = pc[5]; p7 = pc[6]; p8 = pc[7]; p9 = pc[8]; |
77 | p10 = pc[9]; p11 = pc[10]; p12 = pc[11]; p13 = pc[12]; p14 = pc[13]; |
78 | p15 = pc[14]; p16 = pc[15]; p17 = pc[16]; p18 = pc[17]; p19 = pc[18]; |
79 | p20 = pc[19]; p21 = pc[20]; p22 = pc[21]; p23 = pc[22]; p24 = pc[23]; |
80 | p25 = pc[24]; |
81 | if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0 || p5 != 0.0 || |
82 | p6 != 0.0 || p7 != 0.0 || p8 != 0.0 || p9 != 0.0 || p10 != 0.0 || |
83 | p11 != 0.0 || p12 != 0.0 || p13 != 0.0 || p14 != 0.0 || p15 != 0.0 |
84 | || p16 != 0.0 || p17 != 0.0 || p18 != 0.0 || p19 != 0.0 || |
85 | p20 != 0.0 || p21 != 0.0 || p22 != 0.0 || p23 != 0.0 || |
86 | p24 != 0.0 || p25 != 0.0) { |
87 | pv = ba + 25*diag_offset[row]; |
88 | pj = bj + diag_offset[row] + 1; |
89 | x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; |
90 | x5 = pv[4]; x6 = pv[5]; x7 = pv[6]; x8 = pv[7]; x9 = pv[8]; |
91 | x10 = pv[9]; x11 = pv[10]; x12 = pv[11]; x13 = pv[12]; x14 = pv[13]; |
92 | x15 = pv[14]; x16 = pv[15]; x17 = pv[16]; x18 = pv[17]; |
93 | x19 = pv[18]; x20 = pv[19]; x21 = pv[20]; x22 = pv[21]; |
94 | x23 = pv[22]; x24 = pv[23]; x25 = pv[24]; |
95 | pc[0] = m1 = p1*x1 + p6*x2 + p11*x3 + p16*x4 + p21*x5; |
96 | pc[1] = m2 = p2*x1 + p7*x2 + p12*x3 + p17*x4 + p22*x5; |
97 | pc[2] = m3 = p3*x1 + p8*x2 + p13*x3 + p18*x4 + p23*x5; |
98 | pc[3] = m4 = p4*x1 + p9*x2 + p14*x3 + p19*x4 + p24*x5; |
99 | pc[4] = m5 = p5*x1 + p10*x2 + p15*x3 + p20*x4 + p25*x5; |
100 | |
101 | pc[5] = m6 = p1*x6 + p6*x7 + p11*x8 + p16*x9 + p21*x10; |
102 | pc[6] = m7 = p2*x6 + p7*x7 + p12*x8 + p17*x9 + p22*x10; |
103 | pc[7] = m8 = p3*x6 + p8*x7 + p13*x8 + p18*x9 + p23*x10; |
104 | pc[8] = m9 = p4*x6 + p9*x7 + p14*x8 + p19*x9 + p24*x10; |
105 | pc[9] = m10 = p5*x6 + p10*x7 + p15*x8 + p20*x9 + p25*x10; |
106 | |
107 | pc[10] = m11 = p1*x11 + p6*x12 + p11*x13 + p16*x14 + p21*x15; |
108 | pc[11] = m12 = p2*x11 + p7*x12 + p12*x13 + p17*x14 + p22*x15; |
109 | pc[12] = m13 = p3*x11 + p8*x12 + p13*x13 + p18*x14 + p23*x15; |
110 | pc[13] = m14 = p4*x11 + p9*x12 + p14*x13 + p19*x14 + p24*x15; |
111 | pc[14] = m15 = p5*x11 + p10*x12 + p15*x13 + p20*x14 + p25*x15; |
112 | |
113 | pc[15] = m16 = p1*x16 + p6*x17 + p11*x18 + p16*x19 + p21*x20; |
114 | pc[16] = m17 = p2*x16 + p7*x17 + p12*x18 + p17*x19 + p22*x20; |
115 | pc[17] = m18 = p3*x16 + p8*x17 + p13*x18 + p18*x19 + p23*x20; |
116 | pc[18] = m19 = p4*x16 + p9*x17 + p14*x18 + p19*x19 + p24*x20; |
117 | pc[19] = m20 = p5*x16 + p10*x17 + p15*x18 + p20*x19 + p25*x20; |
118 | |
119 | pc[20] = m21 = p1*x21 + p6*x22 + p11*x23 + p16*x24 + p21*x25; |
120 | pc[21] = m22 = p2*x21 + p7*x22 + p12*x23 + p17*x24 + p22*x25; |
121 | pc[22] = m23 = p3*x21 + p8*x22 + p13*x23 + p18*x24 + p23*x25; |
122 | pc[23] = m24 = p4*x21 + p9*x22 + p14*x23 + p19*x24 + p24*x25; |
123 | pc[24] = m25 = p5*x21 + p10*x22 + p15*x23 + p20*x24 + p25*x25; |
124 | |
125 | nz = bi[row+1] - diag_offset[row] - 1; |
126 | pv += 25; |
127 | for (j=0; j<nz; j++) { |
128 | x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; |
129 | x5 = pv[4]; x6 = pv[5]; x7 = pv[6]; x8 = pv[7]; x9 = pv[8]; |
130 | x10 = pv[9]; x11 = pv[10]; x12 = pv[11]; x13 = pv[12]; |
131 | x14 = pv[13]; x15 = pv[14]; x16 = pv[15]; x17 = pv[16]; |
132 | x18 = pv[17]; x19 = pv[18]; x20 = pv[19]; x21 = pv[20]; |
133 | x22 = pv[21]; x23 = pv[22]; x24 = pv[23]; x25 = pv[24]; |
134 | x = rtmp + 25*pj[j]; |
135 | x[0] -= m1*x1 + m6*x2 + m11*x3 + m16*x4 + m21*x5; |
136 | x[1] -= m2*x1 + m7*x2 + m12*x3 + m17*x4 + m22*x5; |
137 | x[2] -= m3*x1 + m8*x2 + m13*x3 + m18*x4 + m23*x5; |
138 | x[3] -= m4*x1 + m9*x2 + m14*x3 + m19*x4 + m24*x5; |
139 | x[4] -= m5*x1 + m10*x2 + m15*x3 + m20*x4 + m25*x5; |
140 | |
141 | x[5] -= m1*x6 + m6*x7 + m11*x8 + m16*x9 + m21*x10; |
142 | x[6] -= m2*x6 + m7*x7 + m12*x8 + m17*x9 + m22*x10; |
143 | x[7] -= m3*x6 + m8*x7 + m13*x8 + m18*x9 + m23*x10; |
144 | x[8] -= m4*x6 + m9*x7 + m14*x8 + m19*x9 + m24*x10; |
145 | x[9] -= m5*x6 + m10*x7 + m15*x8 + m20*x9 + m25*x10; |
146 | |
147 | x[10] -= m1*x11 + m6*x12 + m11*x13 + m16*x14 + m21*x15; |
148 | x[11] -= m2*x11 + m7*x12 + m12*x13 + m17*x14 + m22*x15; |
149 | x[12] -= m3*x11 + m8*x12 + m13*x13 + m18*x14 + m23*x15; |
150 | x[13] -= m4*x11 + m9*x12 + m14*x13 + m19*x14 + m24*x15; |
151 | x[14] -= m5*x11 + m10*x12 + m15*x13 + m20*x14 + m25*x15; |
152 | |
153 | x[15] -= m1*x16 + m6*x17 + m11*x18 + m16*x19 + m21*x20; |
154 | x[16] -= m2*x16 + m7*x17 + m12*x18 + m17*x19 + m22*x20; |
155 | x[17] -= m3*x16 + m8*x17 + m13*x18 + m18*x19 + m23*x20; |
156 | x[18] -= m4*x16 + m9*x17 + m14*x18 + m19*x19 + m24*x20; |
157 | x[19] -= m5*x16 + m10*x17 + m15*x18 + m20*x19 + m25*x20; |
158 | |
159 | x[20] -= m1*x21 + m6*x22 + m11*x23 + m16*x24 + m21*x25; |
160 | x[21] -= m2*x21 + m7*x22 + m12*x23 + m17*x24 + m22*x25; |
161 | x[22] -= m3*x21 + m8*x22 + m13*x23 + m18*x24 + m23*x25; |
162 | x[23] -= m4*x21 + m9*x22 + m14*x23 + m19*x24 + m24*x25; |
163 | x[24] -= m5*x21 + m10*x22 + m15*x23 + m20*x24 + m25*x25; |
164 | |
165 | pv += 25; |
166 | } |
167 | ierr = PetscLogFlops(250.0*nz+225.0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),167,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
168 | } |
169 | row = *ajtmp++; |
170 | } |
171 | /* finished row so stick it into b->a */ |
172 | pv = ba + 25*bi[i]; |
173 | pj = bj + bi[i]; |
174 | nz = bi[i+1] - bi[i]; |
175 | for (j=0; j<nz; j++) { |
176 | #if defined(PETSC_USE_MEMCPY1) |
177 | ierr = PetscArraycpy(pv,rtmp+25*pj[j],25)((sizeof(*(pv)) != sizeof(*(rtmp+25*pj[j]))) || PetscMemcpy(pv ,rtmp+25*pj[j],(25)*sizeof(*(pv))));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),177,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
178 | #else |
179 | x = rtmp+25*pj[j]; |
180 | pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3]; |
181 | pv[4] = x[4]; pv[5] = x[5]; pv[6] = x[6]; pv[7] = x[7]; pv[8] = x[8]; |
182 | pv[9] = x[9]; pv[10] = x[10]; pv[11] = x[11]; pv[12] = x[12]; |
183 | pv[13] = x[13]; pv[14] = x[14]; pv[15] = x[15]; pv[16] = x[16]; |
184 | pv[17] = x[17]; pv[18] = x[18]; pv[19] = x[19]; pv[20] = x[20]; |
185 | pv[21] = x[21]; pv[22] = x[22]; pv[23] = x[23]; pv[24] = x[24]; |
186 | #endif |
187 | pv += 25; |
188 | } |
189 | /* invert diagonal block */ |
190 | w = ba + 25*diag_offset[i]; |
191 | ierr = PetscKernel_A_gets_inverse_A_5(w,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),191,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
192 | if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; |
193 | } |
194 | |
195 | ierr = PetscFree(rtmp)((*PetscTrFree)((void*)(rtmp),195,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ) || ((rtmp) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),195,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
196 | ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),196,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
197 | ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),197,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
198 | |
199 | C->ops->solve = MatSolve_SeqBAIJ_5_inplace; |
200 | C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_inplace; |
201 | C->assembled = PETSC_TRUE; |
202 | |
203 | ierr = PetscLogFlops(1.333333333333*5*5*5*b->mbs);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),203,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */ |
204 | 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); |
205 | } |
206 | |
207 | /* MatLUFactorNumeric_SeqBAIJ_5 - |
208 | copied from MatLUFactorNumeric_SeqBAIJ_N_inplace() and manually re-implemented |
209 | PetscKernel_A_gets_A_times_B() |
210 | PetscKernel_A_gets_A_minus_B_times_C() |
211 | PetscKernel_A_gets_inverse_A() |
212 | */ |
213 | |
214 | PetscErrorCode MatLUFactorNumeric_SeqBAIJ_5(Mat B,Mat A,const MatFactorInfo *info) |
215 | { |
216 | Mat C =B; |
217 | Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data; |
218 | IS isrow = b->row,isicol = b->icol; |
219 | PetscErrorCode ierr; |
220 | const PetscInt *r,*ic; |
221 | PetscInt i,j,k,nz,nzL,row; |
222 | const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j; |
223 | const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag,*pj,bs2=a->bs2; |
224 | MatScalar *rtmp,*pc,*mwork,*v,*pv,*aa=a->a,work[25]; |
225 | PetscInt flg,ipvt[5]; |
226 | PetscReal shift = info->shiftamount; |
227 | PetscBool allowzeropivot,zeropivotdetected; |
228 | |
229 | 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/baijfact9.c" ; petscstack->line[petscstack->currentsize] = 229; petscstack ->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack ->currentsize++; } if (petscstack) { petscstack->hotdepth += (PETSC_FALSE || petscstack->hotdepth); } ; } while (0) ; ; } while (0); |
230 | allowzeropivot = PetscNot(A->erroriffailure)((A->erroriffailure) ? PETSC_FALSE : PETSC_TRUE); |
231 | ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),231,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
232 | ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),232,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
233 | |
234 | /* generate work space needed by the factorization */ |
235 | ierr = PetscMalloc2(bs2*n,&rtmp,bs2,&mwork)PetscMallocA(2,PETSC_FALSE,235,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.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),235,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
236 | ierr = PetscArrayzero(rtmp,bs2*n)PetscMemzero(rtmp,(bs2*n)*sizeof(*(rtmp)));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),236,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
237 | |
238 | for (i=0; i<n; i++) { |
239 | /* zero rtmp */ |
240 | /* L part */ |
241 | nz = bi[i+1] - bi[i]; |
242 | bjtmp = bj + bi[i]; |
243 | for (j=0; j<nz; j++) { |
244 | 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),244,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
245 | } |
246 | |
247 | /* U part */ |
248 | nz = bdiag[i] - bdiag[i+1]; |
249 | bjtmp = bj + bdiag[i+1]+1; |
250 | for (j=0; j<nz; j++) { |
251 | 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),251,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
252 | } |
253 | |
254 | /* load in initial (unfactored row) */ |
255 | nz = ai[r[i]+1] - ai[r[i]]; |
256 | ajtmp = aj + ai[r[i]]; |
257 | v = aa + bs2*ai[r[i]]; |
258 | for (j=0; j<nz; j++) { |
259 | 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),259,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
260 | } |
261 | |
262 | /* elimination */ |
263 | bjtmp = bj + bi[i]; |
264 | nzL = bi[i+1] - bi[i]; |
265 | for (k=0; k < nzL; k++) { |
266 | row = bjtmp[k]; |
267 | pc = rtmp + bs2*row; |
268 | for (flg=0,j=0; j<bs2; j++) { |
269 | if (pc[j]!=0.0) { |
270 | flg = 1; |
271 | break; |
272 | } |
273 | } |
274 | if (flg) { |
275 | pv = b->a + bs2*bdiag[row]; |
276 | /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ |
277 | ierr = PetscKernel_A_gets_A_times_B_5(pc,pv,mwork)0; { ierr = ((sizeof(*(mwork)) != sizeof(*(pc))) || PetscMemcpy (mwork,pc,(25)*sizeof(*(mwork))));;do {if (__builtin_expect(! !(ierr),0)) return PetscError(((MPI_Comm)0x44000001),277,__func__ ,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); pc[0] = mwork[0]*pv [0] + mwork[5]*pv[1] + mwork[10]*pv[2] + mwork[15]*pv[3] + mwork [20]*pv[4]; pc[1] = mwork[1]*pv[0] + mwork[6]*pv[1] + mwork[11 ]*pv[2] + mwork[16]*pv[3] + mwork[21]*pv[4]; pc[2] = mwork[2] *pv[0] + mwork[7]*pv[1] + mwork[12]*pv[2] + mwork[17]*pv[3] + mwork[22]*pv[4]; pc[3] = mwork[3]*pv[0] + mwork[8]*pv[1] + mwork [13]*pv[2] + mwork[18]*pv[3] + mwork[23]*pv[4]; pc[4] = mwork [4]*pv[0] + mwork[9]*pv[1] + mwork[14]*pv[2] + mwork[19]*pv[3 ] + mwork[24]*pv[4]; pc[5] = mwork[0]*pv[5] + mwork[5]*pv[6] + mwork[10]*pv[7] + mwork[15]*pv[8] + mwork[20]*pv[9]; pc[6] = mwork[1]*pv[5] + mwork[6]*pv[6] + mwork[11]*pv[7] + mwork[16 ]*pv[8] + mwork[21]*pv[9]; pc[7] = mwork[2]*pv[5] + mwork[7]* pv[6] + mwork[12]*pv[7] + mwork[17]*pv[8] + mwork[22]*pv[9]; pc [8] = mwork[3]*pv[5] + mwork[8]*pv[6] + mwork[13]*pv[7] + mwork [18]*pv[8] + mwork[23]*pv[9]; pc[9] = mwork[4]*pv[5] + mwork[ 9]*pv[6] + mwork[14]*pv[7] + mwork[19]*pv[8] + mwork[24]*pv[9 ]; pc[10] = mwork[0]*pv[10] + mwork[5]*pv[11] + mwork[10]*pv[ 12] + mwork[15]*pv[13] + mwork[20]*pv[14]; pc[11] = mwork[1]* pv[10] + mwork[6]*pv[11] + mwork[11]*pv[12] + mwork[16]*pv[13 ] + mwork[21]*pv[14]; pc[12] = mwork[2]*pv[10] + mwork[7]*pv[ 11] + mwork[12]*pv[12] + mwork[17]*pv[13] + mwork[22]*pv[14]; pc[13] = mwork[3]*pv[10] + mwork[8]*pv[11] + mwork[13]*pv[12 ] + mwork[18]*pv[13] + mwork[23]*pv[14]; pc[14] = mwork[4]*pv [10] + mwork[9]*pv[11] + mwork[14]*pv[12] + mwork[19]*pv[13] + mwork[24]*pv[14]; pc[15] = mwork[0]*pv[15] + mwork[5]*pv[16] + mwork[10]*pv[17] + mwork[15]*pv[18] + mwork[20]*pv[19]; pc [16] = mwork[1]*pv[15] + mwork[6]*pv[16] + mwork[11]*pv[17] + mwork[16]*pv[18] + mwork[21]*pv[19]; pc[17] = mwork[2]*pv[15 ] + mwork[7]*pv[16] + mwork[12]*pv[17] + mwork[17]*pv[18] + mwork [22]*pv[19]; pc[18] = mwork[3]*pv[15] + mwork[8]*pv[16] + mwork [13]*pv[17] + mwork[18]*pv[18] + mwork[23]*pv[19]; pc[19] = mwork [4]*pv[15] + mwork[9]*pv[16] + mwork[14]*pv[17] + mwork[19]*pv [18] + mwork[24]*pv[19]; pc[20] = mwork[0]*pv[20] + mwork[5]* pv[21] + mwork[10]*pv[22] + mwork[15]*pv[23] + mwork[20]*pv[24 ]; pc[21] = mwork[1]*pv[20] + mwork[6]*pv[21] + mwork[11]*pv[ 22] + mwork[16]*pv[23] + mwork[21]*pv[24]; pc[22] = mwork[2]* pv[20] + mwork[7]*pv[21] + mwork[12]*pv[22] + mwork[17]*pv[23 ] + mwork[22]*pv[24]; pc[23] = mwork[3]*pv[20] + mwork[8]*pv[ 21] + mwork[13]*pv[22] + mwork[18]*pv[23] + mwork[23]*pv[24]; pc[24] = mwork[4]*pv[20] + mwork[9]*pv[21] + mwork[14]*pv[22 ] + mwork[19]*pv[23] + mwork[24]*pv[24]; };CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),277,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
Value stored to 'ierr' is never read | |
278 | |
279 | pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */ |
280 | pv = b->a + bs2*(bdiag[row+1]+1); |
281 | nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */ |
282 | for (j=0; j<nz; j++) { |
283 | /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ |
284 | /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ |
285 | v = rtmp + bs2*pj[j]; |
286 | ierr = PetscKernel_A_gets_A_minus_B_times_C_5(v,pc,pv)0; { v[0] -= pc[0]*pv[0] + pc[5]*pv[1] + pc[10]*pv[2] + pc[15 ]*pv[3] + pc[20]*pv[4]; v[1] -= pc[1]*pv[0] + pc[6]*pv[1] + pc [11]*pv[2] + pc[16]*pv[3] + pc[21]*pv[4]; v[2] -= pc[2]*pv[0] + pc[7]*pv[1] + pc[12]*pv[2] + pc[17]*pv[3] + pc[22]*pv[4]; v [3] -= pc[3]*pv[0] + pc[8]*pv[1] + pc[13]*pv[2] + pc[18]*pv[3 ] + pc[23]*pv[4]; v[4] -= pc[4]*pv[0] + pc[9]*pv[1] + pc[14]* pv[2] + pc[19]*pv[3] + pc[24]*pv[4]; v[5] -= pc[0]*pv[5] + pc [5]*pv[6] + pc[10]*pv[7] + pc[15]*pv[8] + pc[20]*pv[9]; v[6] -= pc[1]*pv[5] + pc[6]*pv[6] + pc[11]*pv[7] + pc[16]*pv[8] + pc [21]*pv[9]; v[7] -= pc[2]*pv[5] + pc[7]*pv[6] + pc[12]*pv[7] + pc[17]*pv[8] + pc[22]*pv[9]; v[8] -= pc[3]*pv[5] + pc[8]*pv[ 6] + pc[13]*pv[7] + pc[18]*pv[8] + pc[23]*pv[9]; v[9] -= pc[4 ]*pv[5] + pc[9]*pv[6] + pc[14]*pv[7] + pc[19]*pv[8] + pc[24]* pv[9]; v[10] -= pc[0]*pv[10] + pc[5]*pv[11] + pc[10]*pv[12] + pc[15]*pv[13] + pc[20]*pv[14]; v[11] -= pc[1]*pv[10] + pc[6] *pv[11] + pc[11]*pv[12] + pc[16]*pv[13] + pc[21]*pv[14]; v[12 ] -= pc[2]*pv[10] + pc[7]*pv[11] + pc[12]*pv[12] + pc[17]*pv[ 13] + pc[22]*pv[14]; v[13] -= pc[3]*pv[10] + pc[8]*pv[11] + pc [13]*pv[12] + pc[18]*pv[13] + pc[23]*pv[14]; v[14] -= pc[4]*pv [10] + pc[9]*pv[11] + pc[14]*pv[12] + pc[19]*pv[13] + pc[24]* pv[14]; v[15] -= pc[0]*pv[15] + pc[5]*pv[16] + pc[10]*pv[17] + pc[15]*pv[18] + pc[20]*pv[19]; v[16] -= pc[1]*pv[15] + pc[6] *pv[16] + pc[11]*pv[17] + pc[16]*pv[18] + pc[21]*pv[19]; v[17 ] -= pc[2]*pv[15] + pc[7]*pv[16] + pc[12]*pv[17] + pc[17]*pv[ 18] + pc[22]*pv[19]; v[18] -= pc[3]*pv[15] + pc[8]*pv[16] + pc [13]*pv[17] + pc[18]*pv[18] + pc[23]*pv[19]; v[19] -= pc[4]*pv [15] + pc[9]*pv[16] + pc[14]*pv[17] + pc[19]*pv[18] + pc[24]* pv[19]; v[20] -= pc[0]*pv[20] + pc[5]*pv[21] + pc[10]*pv[22] + pc[15]*pv[23] + pc[20]*pv[24]; v[21] -= pc[1]*pv[20] + pc[6] *pv[21] + pc[11]*pv[22] + pc[16]*pv[23] + pc[21]*pv[24]; v[22 ] -= pc[2]*pv[20] + pc[7]*pv[21] + pc[12]*pv[22] + pc[17]*pv[ 23] + pc[22]*pv[24]; v[23] -= pc[3]*pv[20] + pc[8]*pv[21] + pc [13]*pv[22] + pc[18]*pv[23] + pc[23]*pv[24]; v[24] -= pc[4]*pv [20] + pc[9]*pv[21] + pc[14]*pv[22] + pc[19]*pv[23] + pc[24]* pv[24]; };CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),286,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
287 | pv += bs2; |
288 | } |
289 | ierr = PetscLogFlops(250*nz+225);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),289,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ |
290 | } |
291 | } |
292 | |
293 | /* finished row so stick it into b->a */ |
294 | /* L part */ |
295 | pv = b->a + bs2*bi[i]; |
296 | pj = b->j + bi[i]; |
297 | nz = bi[i+1] - bi[i]; |
298 | for (j=0; j<nz; j++) { |
299 | 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),299,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
300 | } |
301 | |
302 | /* Mark diagonal and invert diagonal for simplier triangular solves */ |
303 | pv = b->a + bs2*bdiag[i]; |
304 | pj = b->j + bdiag[i]; |
305 | 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),305,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
306 | ierr = PetscKernel_A_gets_inverse_A_5(pv,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),306,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
307 | if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; |
308 | |
309 | /* U part */ |
310 | pv = b->a + bs2*(bdiag[i+1]+1); |
311 | pj = b->j + bdiag[i+1]+1; |
312 | nz = bdiag[i] - bdiag[i+1] - 1; |
313 | for (j=0; j<nz; j++) { |
314 | 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),314,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
315 | } |
316 | } |
317 | |
318 | ierr = PetscFree2(rtmp,mwork)PetscFreeA(2,318,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,&(rtmp),&(mwork));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),318,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
319 | ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),319,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
320 | ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),320,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
321 | |
322 | C->ops->solve = MatSolve_SeqBAIJ_5; |
323 | C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5; |
324 | C->assembled = PETSC_TRUE; |
325 | |
326 | ierr = PetscLogFlops(1.333333333333*5*5*5*n);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),326,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */ |
327 | 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); |
328 | } |
329 | |
330 | /* |
331 | Version for when blocks are 5 by 5 Using natural ordering |
332 | */ |
333 | PetscErrorCode MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info) |
334 | { |
335 | Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data; |
336 | PetscErrorCode ierr; |
337 | PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j,ipvt[5]; |
338 | PetscInt *ajtmpold,*ajtmp,nz,row; |
339 | PetscInt *diag_offset = b->diag,*ai=a->i,*aj=a->j,*pj; |
340 | MatScalar *pv,*v,*rtmp,*pc,*w,*x; |
341 | MatScalar x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15; |
342 | MatScalar x16,x17,x18,x19,x20,x21,x22,x23,x24,x25; |
343 | MatScalar p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15; |
344 | MatScalar p16,p17,p18,p19,p20,p21,p22,p23,p24,p25; |
345 | MatScalar m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15; |
346 | MatScalar m16,m17,m18,m19,m20,m21,m22,m23,m24,m25; |
347 | MatScalar *ba = b->a,*aa = a->a,work[25]; |
348 | PetscReal shift = info->shiftamount; |
349 | PetscBool allowzeropivot,zeropivotdetected; |
350 | |
351 | 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/baijfact9.c" ; petscstack->line[petscstack->currentsize] = 351; petscstack ->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack ->currentsize++; } if (petscstack) { petscstack->hotdepth += (PETSC_FALSE || petscstack->hotdepth); } ; } while (0) ; ; } while (0); |
352 | allowzeropivot = PetscNot(A->erroriffailure)((A->erroriffailure) ? PETSC_FALSE : PETSC_TRUE); |
353 | ierr = PetscMalloc1(25*(n+1),&rtmp)PetscMallocA(1,PETSC_FALSE,353,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,(size_t)(25*(n+1))*sizeof(**(&rtmp)),(&rtmp));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),353,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
354 | for (i=0; i<n; i++) { |
355 | nz = bi[i+1] - bi[i]; |
356 | ajtmp = bj + bi[i]; |
357 | for (j=0; j<nz; j++) { |
358 | x = rtmp+25*ajtmp[j]; |
359 | x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = x[9] = 0.0; |
360 | x[10] = x[11] = x[12] = x[13] = x[14] = x[15] = 0.0; |
361 | x[16] = x[17] = x[18] = x[19] = x[20] = x[21] = x[22] = x[23] = x[24] = 0.0; |
362 | } |
363 | /* load in initial (unfactored row) */ |
364 | nz = ai[i+1] - ai[i]; |
365 | ajtmpold = aj + ai[i]; |
366 | v = aa + 25*ai[i]; |
367 | for (j=0; j<nz; j++) { |
368 | x = rtmp+25*ajtmpold[j]; |
369 | x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3]; |
370 | x[4] = v[4]; x[5] = v[5]; x[6] = v[6]; x[7] = v[7]; x[8] = v[8]; |
371 | x[9] = v[9]; x[10] = v[10]; x[11] = v[11]; x[12] = v[12]; x[13] = v[13]; |
372 | x[14] = v[14]; x[15] = v[15]; x[16] = v[16]; x[17] = v[17]; x[18] = v[18]; |
373 | x[19] = v[19]; x[20] = v[20]; x[21] = v[21]; x[22] = v[22]; x[23] = v[23]; |
374 | x[24] = v[24]; |
375 | v += 25; |
376 | } |
377 | row = *ajtmp++; |
378 | while (row < i) { |
379 | pc = rtmp + 25*row; |
380 | p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3]; |
381 | p5 = pc[4]; p6 = pc[5]; p7 = pc[6]; p8 = pc[7]; p9 = pc[8]; |
382 | p10 = pc[9]; p11 = pc[10]; p12 = pc[11]; p13 = pc[12]; p14 = pc[13]; |
383 | p15 = pc[14]; p16 = pc[15]; p17 = pc[16]; p18 = pc[17]; |
384 | p19 = pc[18]; p20 = pc[19]; p21 = pc[20]; p22 = pc[21]; p23 = pc[22]; |
385 | p24 = pc[23]; p25 = pc[24]; |
386 | if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0 || p5 != 0.0 || |
387 | p6 != 0.0 || p7 != 0.0 || p8 != 0.0 || p9 != 0.0 || p10 != 0.0 || |
388 | p11 != 0.0 || p12 != 0.0 || p13 != 0.0 || p14 != 0.0 || p15 != 0.0 |
389 | || p16 != 0.0 || p17 != 0.0 || p18 != 0.0 || p19 != 0.0 || p20 != 0.0 |
390 | || p21 != 0.0 || p22 != 0.0 || p23 != 0.0 || p24 != 0.0 || p25 != 0.0) { |
391 | pv = ba + 25*diag_offset[row]; |
392 | pj = bj + diag_offset[row] + 1; |
393 | x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; |
394 | x5 = pv[4]; x6 = pv[5]; x7 = pv[6]; x8 = pv[7]; x9 = pv[8]; |
395 | x10 = pv[9]; x11 = pv[10]; x12 = pv[11]; x13 = pv[12]; x14 = pv[13]; |
396 | x15 = pv[14]; x16 = pv[15]; x17 = pv[16]; x18 = pv[17]; x19 = pv[18]; |
397 | x20 = pv[19]; x21 = pv[20]; x22 = pv[21]; x23 = pv[22]; x24 = pv[23]; |
398 | x25 = pv[24]; |
399 | pc[0] = m1 = p1*x1 + p6*x2 + p11*x3 + p16*x4 + p21*x5; |
400 | pc[1] = m2 = p2*x1 + p7*x2 + p12*x3 + p17*x4 + p22*x5; |
401 | pc[2] = m3 = p3*x1 + p8*x2 + p13*x3 + p18*x4 + p23*x5; |
402 | pc[3] = m4 = p4*x1 + p9*x2 + p14*x3 + p19*x4 + p24*x5; |
403 | pc[4] = m5 = p5*x1 + p10*x2 + p15*x3 + p20*x4 + p25*x5; |
404 | |
405 | pc[5] = m6 = p1*x6 + p6*x7 + p11*x8 + p16*x9 + p21*x10; |
406 | pc[6] = m7 = p2*x6 + p7*x7 + p12*x8 + p17*x9 + p22*x10; |
407 | pc[7] = m8 = p3*x6 + p8*x7 + p13*x8 + p18*x9 + p23*x10; |
408 | pc[8] = m9 = p4*x6 + p9*x7 + p14*x8 + p19*x9 + p24*x10; |
409 | pc[9] = m10 = p5*x6 + p10*x7 + p15*x8 + p20*x9 + p25*x10; |
410 | |
411 | pc[10] = m11 = p1*x11 + p6*x12 + p11*x13 + p16*x14 + p21*x15; |
412 | pc[11] = m12 = p2*x11 + p7*x12 + p12*x13 + p17*x14 + p22*x15; |
413 | pc[12] = m13 = p3*x11 + p8*x12 + p13*x13 + p18*x14 + p23*x15; |
414 | pc[13] = m14 = p4*x11 + p9*x12 + p14*x13 + p19*x14 + p24*x15; |
415 | pc[14] = m15 = p5*x11 + p10*x12 + p15*x13 + p20*x14 + p25*x15; |
416 | |
417 | pc[15] = m16 = p1*x16 + p6*x17 + p11*x18 + p16*x19 + p21*x20; |
418 | pc[16] = m17 = p2*x16 + p7*x17 + p12*x18 + p17*x19 + p22*x20; |
419 | pc[17] = m18 = p3*x16 + p8*x17 + p13*x18 + p18*x19 + p23*x20; |
420 | pc[18] = m19 = p4*x16 + p9*x17 + p14*x18 + p19*x19 + p24*x20; |
421 | pc[19] = m20 = p5*x16 + p10*x17 + p15*x18 + p20*x19 + p25*x20; |
422 | |
423 | pc[20] = m21 = p1*x21 + p6*x22 + p11*x23 + p16*x24 + p21*x25; |
424 | pc[21] = m22 = p2*x21 + p7*x22 + p12*x23 + p17*x24 + p22*x25; |
425 | pc[22] = m23 = p3*x21 + p8*x22 + p13*x23 + p18*x24 + p23*x25; |
426 | pc[23] = m24 = p4*x21 + p9*x22 + p14*x23 + p19*x24 + p24*x25; |
427 | pc[24] = m25 = p5*x21 + p10*x22 + p15*x23 + p20*x24 + p25*x25; |
428 | |
429 | nz = bi[row+1] - diag_offset[row] - 1; |
430 | pv += 25; |
431 | for (j=0; j<nz; j++) { |
432 | x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; |
433 | x5 = pv[4]; x6 = pv[5]; x7 = pv[6]; x8 = pv[7]; x9 = pv[8]; |
434 | x10 = pv[9]; x11 = pv[10]; x12 = pv[11]; x13 = pv[12]; |
435 | x14 = pv[13]; x15 = pv[14]; x16 = pv[15]; x17 = pv[16]; x18 = pv[17]; |
436 | x19 = pv[18]; x20 = pv[19]; x21 = pv[20]; x22 = pv[21]; x23 = pv[22]; |
437 | x24 = pv[23]; x25 = pv[24]; |
438 | x = rtmp + 25*pj[j]; |
439 | x[0] -= m1*x1 + m6*x2 + m11*x3 + m16*x4 + m21*x5; |
440 | x[1] -= m2*x1 + m7*x2 + m12*x3 + m17*x4 + m22*x5; |
441 | x[2] -= m3*x1 + m8*x2 + m13*x3 + m18*x4 + m23*x5; |
442 | x[3] -= m4*x1 + m9*x2 + m14*x3 + m19*x4 + m24*x5; |
443 | x[4] -= m5*x1 + m10*x2 + m15*x3 + m20*x4 + m25*x5; |
444 | |
445 | x[5] -= m1*x6 + m6*x7 + m11*x8 + m16*x9 + m21*x10; |
446 | x[6] -= m2*x6 + m7*x7 + m12*x8 + m17*x9 + m22*x10; |
447 | x[7] -= m3*x6 + m8*x7 + m13*x8 + m18*x9 + m23*x10; |
448 | x[8] -= m4*x6 + m9*x7 + m14*x8 + m19*x9 + m24*x10; |
449 | x[9] -= m5*x6 + m10*x7 + m15*x8 + m20*x9 + m25*x10; |
450 | |
451 | x[10] -= m1*x11 + m6*x12 + m11*x13 + m16*x14 + m21*x15; |
452 | x[11] -= m2*x11 + m7*x12 + m12*x13 + m17*x14 + m22*x15; |
453 | x[12] -= m3*x11 + m8*x12 + m13*x13 + m18*x14 + m23*x15; |
454 | x[13] -= m4*x11 + m9*x12 + m14*x13 + m19*x14 + m24*x15; |
455 | x[14] -= m5*x11 + m10*x12 + m15*x13 + m20*x14 + m25*x15; |
456 | |
457 | x[15] -= m1*x16 + m6*x17 + m11*x18 + m16*x19 + m21*x20; |
458 | x[16] -= m2*x16 + m7*x17 + m12*x18 + m17*x19 + m22*x20; |
459 | x[17] -= m3*x16 + m8*x17 + m13*x18 + m18*x19 + m23*x20; |
460 | x[18] -= m4*x16 + m9*x17 + m14*x18 + m19*x19 + m24*x20; |
461 | x[19] -= m5*x16 + m10*x17 + m15*x18 + m20*x19 + m25*x20; |
462 | |
463 | x[20] -= m1*x21 + m6*x22 + m11*x23 + m16*x24 + m21*x25; |
464 | x[21] -= m2*x21 + m7*x22 + m12*x23 + m17*x24 + m22*x25; |
465 | x[22] -= m3*x21 + m8*x22 + m13*x23 + m18*x24 + m23*x25; |
466 | x[23] -= m4*x21 + m9*x22 + m14*x23 + m19*x24 + m24*x25; |
467 | x[24] -= m5*x21 + m10*x22 + m15*x23 + m20*x24 + m25*x25; |
468 | pv += 25; |
469 | } |
470 | ierr = PetscLogFlops(250.0*nz+225.0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),470,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
471 | } |
472 | row = *ajtmp++; |
473 | } |
474 | /* finished row so stick it into b->a */ |
475 | pv = ba + 25*bi[i]; |
476 | pj = bj + bi[i]; |
477 | nz = bi[i+1] - bi[i]; |
478 | for (j=0; j<nz; j++) { |
479 | x = rtmp+25*pj[j]; |
480 | pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3]; |
481 | pv[4] = x[4]; pv[5] = x[5]; pv[6] = x[6]; pv[7] = x[7]; pv[8] = x[8]; |
482 | pv[9] = x[9]; pv[10] = x[10]; pv[11] = x[11]; pv[12] = x[12]; |
483 | pv[13] = x[13]; pv[14] = x[14]; pv[15] = x[15]; pv[16] = x[16]; pv[17] = x[17]; |
484 | pv[18] = x[18]; pv[19] = x[19]; pv[20] = x[20]; pv[21] = x[21]; pv[22] = x[22]; |
485 | pv[23] = x[23]; pv[24] = x[24]; |
486 | pv += 25; |
487 | } |
488 | /* invert diagonal block */ |
489 | w = ba + 25*diag_offset[i]; |
490 | ierr = PetscKernel_A_gets_inverse_A_5(w,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),490,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
491 | if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; |
492 | } |
493 | |
494 | ierr = PetscFree(rtmp)((*PetscTrFree)((void*)(rtmp),494,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ) || ((rtmp) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),494,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
495 | |
496 | C->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering_inplace; |
497 | C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering_inplace; |
498 | C->assembled = PETSC_TRUE; |
499 | |
500 | ierr = PetscLogFlops(1.333333333333*5*5*5*b->mbs);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),500,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */ |
501 | 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); |
502 | } |
503 | |
504 | PetscErrorCode MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info) |
505 | { |
506 | Mat C =B; |
507 | Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data; |
508 | PetscErrorCode ierr; |
509 | PetscInt i,j,k,nz,nzL,row; |
510 | const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j; |
511 | const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag,*pj,bs2=a->bs2; |
512 | MatScalar *rtmp,*pc,*mwork,*v,*vv,*pv,*aa=a->a,work[25]; |
513 | PetscInt flg,ipvt[5]; |
514 | PetscReal shift = info->shiftamount; |
515 | PetscBool allowzeropivot,zeropivotdetected; |
516 | |
517 | 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/baijfact9.c" ; petscstack->line[petscstack->currentsize] = 517; petscstack ->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack ->currentsize++; } if (petscstack) { petscstack->hotdepth += (PETSC_FALSE || petscstack->hotdepth); } ; } while (0) ; ; } while (0); |
518 | allowzeropivot = PetscNot(A->erroriffailure)((A->erroriffailure) ? PETSC_FALSE : PETSC_TRUE); |
519 | |
520 | /* generate work space needed by the factorization */ |
521 | ierr = PetscMalloc2(bs2*n,&rtmp,bs2,&mwork)PetscMallocA(2,PETSC_FALSE,521,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.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),521,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
522 | ierr = PetscArrayzero(rtmp,bs2*n)PetscMemzero(rtmp,(bs2*n)*sizeof(*(rtmp)));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),522,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
523 | |
524 | for (i=0; i<n; i++) { |
525 | /* zero rtmp */ |
526 | /* L part */ |
527 | nz = bi[i+1] - bi[i]; |
528 | bjtmp = bj + bi[i]; |
529 | for (j=0; j<nz; j++) { |
530 | 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),530,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
531 | } |
532 | |
533 | /* U part */ |
534 | nz = bdiag[i] - bdiag[i+1]; |
535 | bjtmp = bj + bdiag[i+1]+1; |
536 | for (j=0; j<nz; j++) { |
537 | 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),537,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
538 | } |
539 | |
540 | /* load in initial (unfactored row) */ |
541 | nz = ai[i+1] - ai[i]; |
542 | ajtmp = aj + ai[i]; |
543 | v = aa + bs2*ai[i]; |
544 | for (j=0; j<nz; j++) { |
545 | 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),545,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
546 | } |
547 | |
548 | /* elimination */ |
549 | bjtmp = bj + bi[i]; |
550 | nzL = bi[i+1] - bi[i]; |
551 | for (k=0; k < nzL; k++) { |
552 | row = bjtmp[k]; |
553 | pc = rtmp + bs2*row; |
554 | for (flg=0,j=0; j<bs2; j++) { |
555 | if (pc[j]!=0.0) { |
556 | flg = 1; |
557 | break; |
558 | } |
559 | } |
560 | if (flg) { |
561 | pv = b->a + bs2*bdiag[row]; |
562 | /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ |
563 | ierr = PetscKernel_A_gets_A_times_B_5(pc,pv,mwork)0; { ierr = ((sizeof(*(mwork)) != sizeof(*(pc))) || PetscMemcpy (mwork,pc,(25)*sizeof(*(mwork))));;do {if (__builtin_expect(! !(ierr),0)) return PetscError(((MPI_Comm)0x44000001),563,__func__ ,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); pc[0] = mwork[0]*pv [0] + mwork[5]*pv[1] + mwork[10]*pv[2] + mwork[15]*pv[3] + mwork [20]*pv[4]; pc[1] = mwork[1]*pv[0] + mwork[6]*pv[1] + mwork[11 ]*pv[2] + mwork[16]*pv[3] + mwork[21]*pv[4]; pc[2] = mwork[2] *pv[0] + mwork[7]*pv[1] + mwork[12]*pv[2] + mwork[17]*pv[3] + mwork[22]*pv[4]; pc[3] = mwork[3]*pv[0] + mwork[8]*pv[1] + mwork [13]*pv[2] + mwork[18]*pv[3] + mwork[23]*pv[4]; pc[4] = mwork [4]*pv[0] + mwork[9]*pv[1] + mwork[14]*pv[2] + mwork[19]*pv[3 ] + mwork[24]*pv[4]; pc[5] = mwork[0]*pv[5] + mwork[5]*pv[6] + mwork[10]*pv[7] + mwork[15]*pv[8] + mwork[20]*pv[9]; pc[6] = mwork[1]*pv[5] + mwork[6]*pv[6] + mwork[11]*pv[7] + mwork[16 ]*pv[8] + mwork[21]*pv[9]; pc[7] = mwork[2]*pv[5] + mwork[7]* pv[6] + mwork[12]*pv[7] + mwork[17]*pv[8] + mwork[22]*pv[9]; pc [8] = mwork[3]*pv[5] + mwork[8]*pv[6] + mwork[13]*pv[7] + mwork [18]*pv[8] + mwork[23]*pv[9]; pc[9] = mwork[4]*pv[5] + mwork[ 9]*pv[6] + mwork[14]*pv[7] + mwork[19]*pv[8] + mwork[24]*pv[9 ]; pc[10] = mwork[0]*pv[10] + mwork[5]*pv[11] + mwork[10]*pv[ 12] + mwork[15]*pv[13] + mwork[20]*pv[14]; pc[11] = mwork[1]* pv[10] + mwork[6]*pv[11] + mwork[11]*pv[12] + mwork[16]*pv[13 ] + mwork[21]*pv[14]; pc[12] = mwork[2]*pv[10] + mwork[7]*pv[ 11] + mwork[12]*pv[12] + mwork[17]*pv[13] + mwork[22]*pv[14]; pc[13] = mwork[3]*pv[10] + mwork[8]*pv[11] + mwork[13]*pv[12 ] + mwork[18]*pv[13] + mwork[23]*pv[14]; pc[14] = mwork[4]*pv [10] + mwork[9]*pv[11] + mwork[14]*pv[12] + mwork[19]*pv[13] + mwork[24]*pv[14]; pc[15] = mwork[0]*pv[15] + mwork[5]*pv[16] + mwork[10]*pv[17] + mwork[15]*pv[18] + mwork[20]*pv[19]; pc [16] = mwork[1]*pv[15] + mwork[6]*pv[16] + mwork[11]*pv[17] + mwork[16]*pv[18] + mwork[21]*pv[19]; pc[17] = mwork[2]*pv[15 ] + mwork[7]*pv[16] + mwork[12]*pv[17] + mwork[17]*pv[18] + mwork [22]*pv[19]; pc[18] = mwork[3]*pv[15] + mwork[8]*pv[16] + mwork [13]*pv[17] + mwork[18]*pv[18] + mwork[23]*pv[19]; pc[19] = mwork [4]*pv[15] + mwork[9]*pv[16] + mwork[14]*pv[17] + mwork[19]*pv [18] + mwork[24]*pv[19]; pc[20] = mwork[0]*pv[20] + mwork[5]* pv[21] + mwork[10]*pv[22] + mwork[15]*pv[23] + mwork[20]*pv[24 ]; pc[21] = mwork[1]*pv[20] + mwork[6]*pv[21] + mwork[11]*pv[ 22] + mwork[16]*pv[23] + mwork[21]*pv[24]; pc[22] = mwork[2]* pv[20] + mwork[7]*pv[21] + mwork[12]*pv[22] + mwork[17]*pv[23 ] + mwork[22]*pv[24]; pc[23] = mwork[3]*pv[20] + mwork[8]*pv[ 21] + mwork[13]*pv[22] + mwork[18]*pv[23] + mwork[23]*pv[24]; pc[24] = mwork[4]*pv[20] + mwork[9]*pv[21] + mwork[14]*pv[22 ] + mwork[19]*pv[23] + mwork[24]*pv[24]; };CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),563,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
564 | |
565 | pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */ |
566 | pv = b->a + bs2*(bdiag[row+1]+1); |
567 | nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */ |
568 | for (j=0; j<nz; j++) { |
569 | /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ |
570 | /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ |
571 | vv = rtmp + bs2*pj[j]; |
572 | ierr = PetscKernel_A_gets_A_minus_B_times_C_5(vv,pc,pv)0; { vv[0] -= pc[0]*pv[0] + pc[5]*pv[1] + pc[10]*pv[2] + pc[15 ]*pv[3] + pc[20]*pv[4]; vv[1] -= pc[1]*pv[0] + pc[6]*pv[1] + pc [11]*pv[2] + pc[16]*pv[3] + pc[21]*pv[4]; vv[2] -= pc[2]*pv[0 ] + pc[7]*pv[1] + pc[12]*pv[2] + pc[17]*pv[3] + pc[22]*pv[4]; vv[3] -= pc[3]*pv[0] + pc[8]*pv[1] + pc[13]*pv[2] + pc[18]*pv [3] + pc[23]*pv[4]; vv[4] -= pc[4]*pv[0] + pc[9]*pv[1] + pc[14 ]*pv[2] + pc[19]*pv[3] + pc[24]*pv[4]; vv[5] -= pc[0]*pv[5] + pc[5]*pv[6] + pc[10]*pv[7] + pc[15]*pv[8] + pc[20]*pv[9]; vv [6] -= pc[1]*pv[5] + pc[6]*pv[6] + pc[11]*pv[7] + pc[16]*pv[8 ] + pc[21]*pv[9]; vv[7] -= pc[2]*pv[5] + pc[7]*pv[6] + pc[12] *pv[7] + pc[17]*pv[8] + pc[22]*pv[9]; vv[8] -= pc[3]*pv[5] + pc [8]*pv[6] + pc[13]*pv[7] + pc[18]*pv[8] + pc[23]*pv[9]; vv[9] -= pc[4]*pv[5] + pc[9]*pv[6] + pc[14]*pv[7] + pc[19]*pv[8] + pc[24]*pv[9]; vv[10] -= pc[0]*pv[10] + pc[5]*pv[11] + pc[10] *pv[12] + pc[15]*pv[13] + pc[20]*pv[14]; vv[11] -= pc[1]*pv[10 ] + pc[6]*pv[11] + pc[11]*pv[12] + pc[16]*pv[13] + pc[21]*pv[ 14]; vv[12] -= pc[2]*pv[10] + pc[7]*pv[11] + pc[12]*pv[12] + pc [17]*pv[13] + pc[22]*pv[14]; vv[13] -= pc[3]*pv[10] + pc[8]*pv [11] + pc[13]*pv[12] + pc[18]*pv[13] + pc[23]*pv[14]; vv[14] -= pc[4]*pv[10] + pc[9]*pv[11] + pc[14]*pv[12] + pc[19]*pv[13] + pc[24]*pv[14]; vv[15] -= pc[0]*pv[15] + pc[5]*pv[16] + pc[10 ]*pv[17] + pc[15]*pv[18] + pc[20]*pv[19]; vv[16] -= pc[1]*pv[ 15] + pc[6]*pv[16] + pc[11]*pv[17] + pc[16]*pv[18] + pc[21]*pv [19]; vv[17] -= pc[2]*pv[15] + pc[7]*pv[16] + pc[12]*pv[17] + pc[17]*pv[18] + pc[22]*pv[19]; vv[18] -= pc[3]*pv[15] + pc[8 ]*pv[16] + pc[13]*pv[17] + pc[18]*pv[18] + pc[23]*pv[19]; vv[ 19] -= pc[4]*pv[15] + pc[9]*pv[16] + pc[14]*pv[17] + pc[19]*pv [18] + pc[24]*pv[19]; vv[20] -= pc[0]*pv[20] + pc[5]*pv[21] + pc[10]*pv[22] + pc[15]*pv[23] + pc[20]*pv[24]; vv[21] -= pc[ 1]*pv[20] + pc[6]*pv[21] + pc[11]*pv[22] + pc[16]*pv[23] + pc [21]*pv[24]; vv[22] -= pc[2]*pv[20] + pc[7]*pv[21] + pc[12]*pv [22] + pc[17]*pv[23] + pc[22]*pv[24]; vv[23] -= pc[3]*pv[20] + pc[8]*pv[21] + pc[13]*pv[22] + pc[18]*pv[23] + pc[23]*pv[24] ; vv[24] -= pc[4]*pv[20] + pc[9]*pv[21] + pc[14]*pv[22] + pc[ 19]*pv[23] + pc[24]*pv[24]; };CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),572,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
573 | pv += bs2; |
574 | } |
575 | ierr = PetscLogFlops(250*nz+225);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),575,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ |
576 | } |
577 | } |
578 | |
579 | /* finished row so stick it into b->a */ |
580 | /* L part */ |
581 | pv = b->a + bs2*bi[i]; |
582 | pj = b->j + bi[i]; |
583 | nz = bi[i+1] - bi[i]; |
584 | for (j=0; j<nz; j++) { |
585 | 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),585,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
586 | } |
587 | |
588 | /* Mark diagonal and invert diagonal for simplier triangular solves */ |
589 | pv = b->a + bs2*bdiag[i]; |
590 | pj = b->j + bdiag[i]; |
591 | 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),591,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
592 | ierr = PetscKernel_A_gets_inverse_A_5(pv,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),592,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
593 | if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; |
594 | |
595 | /* U part */ |
596 | pv = b->a + bs2*(bdiag[i+1]+1); |
597 | pj = b->j + bdiag[i+1]+1; |
598 | nz = bdiag[i] - bdiag[i+1] - 1; |
599 | for (j=0; j<nz; j++) { |
600 | 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),600,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
601 | } |
602 | } |
603 | ierr = PetscFree2(rtmp,mwork)PetscFreeA(2,603,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,&(rtmp),&(mwork));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),603,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); |
604 | |
605 | C->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; |
606 | C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering; |
607 | C->assembled = PETSC_TRUE; |
608 | |
609 | ierr = PetscLogFlops(1.333333333333*5*5*5*n);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),609,__func__,"/sandbox/petsc/petsc.master/src/mat/impls/baij/seq/baijfact9.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */ |
610 | 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); |
611 | } |