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