/sandbox/petsc/petsc.next-tmp/src/mat/impls/baij/seq/baijfact.c

Bug Summary

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

Annotated Source Code

[?] Use j/k keys for keyboard navigation

/sandbox/petsc/petsc.next-tmp/src/mat/impls/baij/seq/baijfact.c

→

Factorization code for BAIJ format.

#include <../src/mat/impls/baij/seq/baij.h>

#include <petsc/private/kernels/blockinvert.h>

PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2(Mat B,Mat A,const MatFactorInfo *info)

{

Mat C =B;

Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data;

IS isrow = b->row,isicol = b->icol;

PetscErrorCode ierr;

const PetscInt *r,*ic;

PetscInt i,j,k,nz,nzL,row,*pj;

const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2;

const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag;

MatScalar *rtmp,*pc,*mwork,*pv;

MatScalar *aa=a->a,*v;

PetscInt flg;

PetscReal shift = info->shiftamount;

PetscBool allowzeropivot,zeropivotdetected;

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/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);

ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),25,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),26,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

allowzeropivot = PetscNot(A->erroriffailure)((A->erroriffailure) ? PETSC_FALSE : PETSC_TRUE);

/* generate work space needed by the factorization */

ierr = PetscMalloc2(bs2*n,&rtmp,bs2,&mwork)PetscMallocA(2,PETSC_FALSE,30,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

for (i=0; i<n; i++) {

/* zero rtmp */

/* L part */

nz = bi[i+1] - bi[i];

bjtmp = bj + bi[i];

for (j=0; j<nz; j++) {

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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

}

/* U part */

nz = bdiag[i] - bdiag[i+1];

bjtmp = bj + bdiag[i+1]+1;

for (j=0; j<nz; j++) {

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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

}

/* load in initial (unfactored row) */

nz = ai[r[i]+1] - ai[r[i]];

ajtmp = aj + ai[r[i]];

v = aa + bs2*ai[r[i]];

for (j=0; j<nz; j++) {

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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

}

/* elimination */

bjtmp = bj + bi[i];

nzL = bi[i+1] - bi[i];

for (k=0; k < nzL; k++) {

row = bjtmp[k];

pc = rtmp + bs2*row;

for (flg=0,j=0; j<bs2; j++) {

if (pc[j] != (PetscScalar)0.0) {

flg = 1;

break;

}

if (flg) {

pv = b->a + bs2*bdiag[row];

/* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */

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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */

pv = b->a + bs2*(bdiag[row+1]+1);

nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */

for (j=0; j<nz; j++) {

/* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */

/* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */

v = rtmp + 4*pj[j];

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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

pv += 4;

}

ierr = PetscLogFlops(16*nz+12);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),84,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */

}

/* finished row so stick it into b->a */

/* L part */

pv = b->a + bs2*bi[i];

pj = b->j + bi[i];

nz = bi[i+1] - bi[i];

for (j=0; j<nz; j++) {

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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

}

/* Mark diagonal and invert diagonal for simplier triangular solves */

pv = b->a + bs2*bdiag[i];

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */

239

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */

331

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* from inverting diagonal blocks */

423

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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

451

452

if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */

453

ddiag = a->diag;

454

sctx.shift_top = info->zeropivot;

455

for (i=0; i<n; i++) {

456

/* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */

457

d = (aa)[ddiag[i]];

458

rs = -PetscAbsScalar(d)fabs(d) - PetscRealPart(d)(d);

459

v = aa+ai[i];

460

nz = ai[i+1] - ai[i];

461

for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j])fabs(v[j]);

462

if (rs>sctx.shift_top) sctx.shift_top = rs;

463

}

464

sctx.shift_top *= 1.1;

465

sctx.nshift_max = 5;

466

sctx.shift_lo = 0.;

467

sctx.shift_hi = 1.;

468

}

469

470

ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),470,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

514

}

515

row = *bjtmp++;

516

}

517

518

/* finished row so stick it into b->a */

519

rs = 0.0;

520

/* L part */

521

pv = b->a + bi[i];

522

pj = b->j + bi[i];

523

nz = bi[i+1] - bi[i];

524

for (j=0; j<nz; j++) {

525

pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j])fabs(pv[j]);

526

}

527

528

/* U part */

529

pv = b->a + bdiag[i+1]+1;

530

pj = b->j + bdiag[i+1]+1;

531

nz = bdiag[i] - bdiag[i+1]-1;

532

for (j=0; j<nz; j++) {

533

pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j])fabs(pv[j]);

534

}

535

536

sctx.rs = rs;

537

sctx.pv = rtmp[i];

538

ierr = MatPivotCheck(B,A,info,&sctx,i);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),538,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

588

}

589

}

590

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

661

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

691

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

710

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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

735

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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

743

744

sctx.shift_amount = 0.;

745

sctx.nshift = 0;

746

do {

747

sctx.newshift = PETSC_FALSE;

748

for (i=0; i<mbs; i++) {

749

rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0;

750

}

751

752

for (k = 0; k<mbs; k++) {

753

bval = ba + bi[k];

754

/* initialize k-th row by the perm[k]-th row of A */

755

jmin = ai[rip[k]]; jmax = ai[rip[k]+1];

756

for (j = jmin; j < jmax; j++) {

757

col = rip[aj[j]];

758

if (col >= k) { /* only take upper triangular entry */

759

rtmp[col] = aa[j];

760

*bval++ = 0.0; /* for in-place factorization */

761

}

762

}

763

764

/* shift the diagonal of the matrix */

765

if (sctx.nshift) rtmp[k] += sctx.shift_amount;

766

767

/* modify k-th row by adding in those rows i with U(i,k)!=0 */

768

dk = rtmp[k];

769

i = jl[k]; /* first row to be added to k_th row */

770

771

while (i < k) {

772

nexti = jl[i]; /* next row to be added to k_th row */

773

774

/* compute multiplier, update diag(k) and U(i,k) */

775

ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */

776

uikdi = -ba[ili]*ba[bi[i]]; /* diagonal(k) */

777

dk += uikdi*ba[ili];

778

ba[ili] = uikdi; /* -U(i,k) */

779

780

/* add multiple of row i to k-th row */

781

jmin = ili + 1; jmax = bi[i+1];

782

if (jmin < jmax) {

783

for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j];

784

/* update il and jl for row i */

785

il[i] = jmin;

786

j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;

787

}

788

i = nexti;

789

}

790

791

/* shift the diagonals when zero pivot is detected */

792

/* compute rs=sum of abs(off-diagonal) */

793

rs = 0.0;

794

jmin = bi[k]+1;

795

nz = bi[k+1] - jmin;

796

if (nz) {

797

bcol = bj + jmin;

798

while (nz--) {

799

rs += PetscAbsScalar(rtmp[*bcol])fabs(rtmp[*bcol]);

800

bcol++;

801

}

802

}

803

804

sctx.rs = rs;

805

sctx.pv = dk;

806

ierr = MatPivotCheck(C,A,info,&sctx,k);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),806,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

836

}

837

}

838

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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

858

859

do {

860

sctx.newshift = PETSC_FALSE;

861

for (i=0; i<am; i++) {

862

rtmp[i] = 0.0; jl[i] = am; il[0] = 0;

863

}

864

865

for (k = 0; k<am; k++) {

866

/* initialize k-th row with elements nonzero in row perm(k) of A */

867

nz = ai[k+1] - ai[k];

868

acol = aj + ai[k];

869

aval = aa + ai[k];

870

bval = ba + bi[k];

871

while (nz--) {

872

if (*acol < k) { /* skip lower triangular entries */

873

acol++; aval++;

874

} else {

875

rtmp[*acol++] = *aval++;

876

*bval++ = 0.0; /* for in-place factorization */

877

}

878

}

879

880

/* shift the diagonal of the matrix */

881

if (sctx.nshift) rtmp[k] += sctx.shift_amount;

882

883

/* modify k-th row by adding in those rows i with U(i,k)!=0 */

884

dk = rtmp[k];

885

i = jl[k]; /* first row to be added to k_th row */

886

887

while (i < k) {

888

nexti = jl[i]; /* next row to be added to k_th row */

889

/* compute multiplier, update D(k) and U(i,k) */

890

ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */

891

uikdi = -ba[ili]*ba[bi[i]];

892

dk += uikdi*ba[ili];

893

ba[ili] = uikdi; /* -U(i,k) */

894

895

/* add multiple of row i to k-th row ... */

896

jmin = ili + 1;

897

nz = bi[i+1] - jmin;

898

if (nz > 0) {

899

bcol = bj + jmin;

900

bval = ba + jmin;

901

while (nz--) rtmp[*bcol++] += uikdi*(*bval++);

902

/* update il and jl for i-th row */

903

il[i] = jmin;

904

j = bj[jmin]; jl[i] = jl[j]; jl[j] = i;

905

}

906

i = nexti;

907

}

908

909

/* shift the diagonals when zero pivot is detected */

910

/* compute rs=sum of abs(off-diagonal) */

911

rs = 0.0;

912

jmin = bi[k]+1;

913

nz = bi[k+1] - jmin;

914

if (nz) {

915

bcol = bj + jmin;

916

while (nz--) {

917

rs += PetscAbsScalar(rtmp[*bcol])fabs(rtmp[*bcol]);

918

bcol++;

919

}

920

}

921

922

sctx.rs = rs;

923

sctx.pv = dk;

924

ierr = MatPivotCheck(C,A,info,&sctx,k);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),924,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

958

}

959

}

960

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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,73,PETSC_ERROR_INITIAL,"Matrix is missing diagonal entry %D"
,i);

Assuming 'missing' is 0

→

←

Taking false branch

→

984

985

if (bs > 1) {

←

Assuming 'bs' is <= 1

→

←

Taking false branch

→

986

if (!a->sbaijMat) {

987

ierr = MatConvert(A,MATSEQSBAIJ"seqsbaij",MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),987,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

992

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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

1023

}

1024

1025

/* initialization */

1026

ierr = PetscMalloc1(am+1,&ui)PetscMallocA(1,PETSC_FALSE,1026,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1033

for (i=0; i<am; i++) {

←

Assuming 'i' is < 'am'

→

←

Loop condition is true. Entering loop body

→

←

Assuming 'i' is >= 'am'

→

←

Loop condition is false. Execution continues on line 1038

→

1034

jl[i] = am; il[i] = 0;

1035

}

1036

1037

/* create and initialize a linked list for storing column indices of the active row k */

1038

nlnk = am + 1;

1039

ierr = PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt)(PetscIntMultError(2,nlnk,((void*)0)) || PetscMallocA(1,PETSC_FALSE
,1039,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/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-tmp/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-tmp/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 */

←

Loop condition is true. Entering loop body

→

1050

/* initialize lnk by the column indices of row rip[k] of A */

1051

nzk = 0;

1052

ncols = ai[rip[k]+1] - ai[rip[k]];

1053

ncols_upper = 0;

1054

cols = cols_lvl + am;

1055

for (j=0; j<ncols; j++) {

←

Assuming 'j' is < 'ncols'

→

←

Loop condition is true. Entering loop body

→

←

Assuming 'j' is < 'ncols'

→

←

Loop condition is true. Entering loop body

→

←

Assuming 'j' is >= 'ncols'

→

←

Loop condition is false. Execution continues on line 1063

→

1056

i = rip[*(aj + ai[rip[k]] + j)];

1057

if (i >= k) { /* only take upper triangular entry */

←

Assuming 'i' is >= 'k'

Taking true branch

Assuming 'i' is >= 'k'

→

←

Taking true branch

→

1058

cols[ncols_upper] = i;

1059

cols_lvl[ncols_upper] = -1; /* initialize level for nonzero entries */

1060

ncols_upper++;

1061

}

1062

}

1063

ierr = PetscIncompleteLLAdd(ncols_upper,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt)0;{ PetscInt _k,_entry,_location,_lnkdata,_incrlev; nlnk = 0;
_lnkdata = am; for (_k=0; _k<ncols_upper; _k++){ _incrlev
= cols_lvl[_k] + 1; if (_incrlev > levels) continue; _entry
= cols[_k]; if (!PetscBTLookupSet(lnkbt,_entry)){ if (_k &&
_entry < _lnkdata) _lnkdata = am; do { _location = _lnkdata
; _lnkdata = lnk[_location]; } while (_entry > _lnkdata); lnk
[_location] = _entry; lnk[_entry] = _lnkdata; lnk_lvl[_entry]
= _incrlev; nlnk++; _lnkdata = _entry; } else { if (lnk_lvl[
_entry] > _incrlev) lnk_lvl[_entry] = _incrlev; } }};CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1063,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

←

Within the expansion of the macro 'PetscIncompleteLLAdd':

→

a	Calling 'PetscBTLookupSet'

1064

nzk += nlnk;

1065

1066

/* update lnk by computing fill-in for each pivot row to be merged in */

1067

prow = jl[k]; /* 1st pivot row */

1068

1069

while (prow < k) {

1070

nextprow = jl[prow];

1071

1072

/* merge prow into k-th row */

1073

jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */

1074

jmax = ui[prow+1];

1075

ncols = jmax-jmin;

1076

i = jmin - ui[prow];

1077

cols = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */

1078

for (j=0; j<ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j);

1079

ierr = PetscIncompleteLLAddSorted(ncols,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt)0;{ PetscInt _k,_entry,_location,_lnkdata,_incrlev; nlnk = 0;
_lnkdata = am; for (_k=0; _k<ncols; _k++){ _incrlev = cols_lvl
[_k] + 1; if (_incrlev > levels) continue; _entry = cols[_k
]; if (!PetscBTLookupSet(lnkbt,_entry)){ do { _location = _lnkdata
; _lnkdata = lnk[_location]; } while (_entry > _lnkdata); lnk
[_location] = _entry; lnk[_entry] = _lnkdata; lnk_lvl[_entry]
= _incrlev; nlnk++; _lnkdata = _entry; } else { if (lnk_lvl[
_entry] > _incrlev) lnk_lvl[_entry] = _incrlev; } }};CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1079,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1080

nzk += nlnk;

1081

1082

/* update il and jl for prow */

1083

if (jmin < jmax) {

1084

il[prow] = jmin;

1085

1086

j = *cols; jl[prow] = jl[j]; jl[j] = prow;

1087

}

1088

prow = nextprow;

1089

}

1090

1091

/* if free space is not available, make more free space */

1092

if (current_space->local_remaining<nzk) {

1093

i = am - k + 1; /* num of unfactored rows */

1094

i = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1))(((PetscIntMultTruncate(i,nzk))<(PetscIntMultTruncate(i,i-
1))) ? (PetscIntMultTruncate(i,nzk)) : (PetscIntMultTruncate(
i,i-1))); /* i*nzk, i*(i-1): estimated and max additional space needed */

1095

ierr = PetscFreeSpaceGet(i,&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1095,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1096

ierr = PetscFreeSpaceGet(i,&current_space_lvl);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1096,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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

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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1213

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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

1286

ierr = PetscFreeSpaceGet(i,&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1286,__func__,"/sandbox/petsc/petsc.next-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1416

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1472

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-tmp/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-tmp/src/mat/impls/baij/seq/baijfact.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1485

for (i=0; i<nbs; i++) {

1486

for (j=0; j<bs2; j++) {

1487

if (absarray[i] < PetscAbsScalar(blockarray[i*nbs+j])fabs(blockarray[i*nbs+j])) absarray[i] = PetscAbsScalar(blockarray[i*nbs+j])fabs(blockarray[i*nbs+j]);

1488

}

1489

}

1490

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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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-tmp/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

1740

}

←

/sandbox/petsc/petsc.next-tmp/include/petscbt.h

#if !defined(PETSCBT_H)

#define PETSCBT_H

#include <petscconf.h>

#include <petscviewer.h>

/*S

PetscBT - PETSc bitarrays

Level: advanced

PetscBTCreate(m,&bt) - creates a bit array with enough room to hold m values

PetscBTDestroy(&bt) - destroys the bit array

PetscBTMemzero(m,bt) - zeros the entire bit array (sets all values to false)

PetscBTSet(bt,index) - sets a particular entry as true

PetscBTClear(bt,index) - sets a particular entry as false

PetscBTLookup(bt,index) - returns the value

PetscBTLookupSet(bt,index) - returns the value and then sets it true

PetscBTLookupClear(bt,index) - returns the value and then sets it false

PetscBTLength(m) - returns number of bytes in array with m bits

PetscBTView(m,bt,viewer) - prints all the entries in a bit array

We do not currently check error flags on PetscBTSet(), PetscBTClear(), PetscBTLookup(),

PetcBTLookupSet(), PetscBTLength() cause error checking would cost hundreds more cycles then

the operation.

S*/

typedef char* PetscBT;

PETSC_STATIC_INLINEstatic inline PetscInt PetscBTLength(PetscInt m)

{

return m/PETSC_BITS_PER_BYTE8+1;

}

PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTMemzero(PetscInt m,PetscBT array)

{

return PetscMemzero(array,sizeof(char)*((size_t)m/PETSC_BITS_PER_BYTE8+1));

}

PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTDestroy(PetscBT *array)

{

return PetscFree(*array)((*PetscTrFree)((void*)(*array),44,__func__,"/sandbox/petsc/petsc.next-tmp/include/petscbt.h"
) || ((*array) = 0,0));

}

PETSC_STATIC_INLINEstatic inline char PetscBTLookup(PetscBT array,PetscInt index)

{

char BT_mask,BT_c;

PetscInt BT_idx;

BT_idx = index/PETSC_BITS_PER_BYTE8;

BT_c = array[BT_idx];

BT_mask = (char)(1 << index%PETSC_BITS_PER_BYTE8);

return (char)(BT_c & BT_mask);

}

PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTView(PetscInt m,const PetscBT bt,PetscViewer viewer)

{

PetscInt i;

PetscErrorCode ierr;

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-tmp/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);}

ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),64,__func__,"/sandbox/petsc/petsc.next-tmp/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

for (i=0; i<m; i++) {

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-tmp/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

}

ierr = PetscViewerFlush(viewer);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),68,__func__,"/sandbox/petsc/petsc.next-tmp/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),69,__func__,"/sandbox/petsc/petsc.next-tmp/include/petscbt.h"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

return 0;

}

PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTCreate(PetscInt m,PetscBT *array)

{

return PetscMalloc1((size_t)m/PETSC_BITS_PER_BYTE+1,array)PetscMallocA(1,PETSC_FALSE,75,__func__,"/sandbox/petsc/petsc.next-tmp/include/petscbt.h"
,(size_t)((size_t)m/8 +1)*sizeof(**(array)),(array)) || PetscBTMemzero(m,*array);

}

PETSC_STATIC_INLINEstatic inline char PetscBTLookupSet(PetscBT array,PetscInt index)

{

char BT_mask,BT_c;

PetscInt BT_idx;

BT_idx = index/PETSC_BITS_PER_BYTE8;

BT_c = array[BT_idx];

BT_mask = (char)(1 << index%PETSC_BITS_PER_BYTE8);

←

The result of the left shift is undefined because the right operand is negative

array[BT_idx] = (char)(BT_c | BT_mask);

return (char)(BT_c & BT_mask);

}

PETSC_STATIC_INLINEstatic inline PetscErrorCode PetscBTSet(PetscBT array,PetscInt index)

{

char BT_mask,BT_c;

PetscInt BT_idx;

BT_idx = index/PETSC_BITS_PER_BYTE8;

BT_c = array[BT_idx];

BT_mask = (char)(1 << index%PETSC_BITS_PER_BYTE8);

array[BT_idx] = (char)(BT_c | BT_mask);

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