/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c

Bug Summary

File:	mat/impls/aij/seq/matmatmult.c
Warning:	line 686, column 33 Dereference of null pointer

Annotated Source Code

[?] Use j/k keys for keyboard navigation

Defines matrix-matrix product routines for pairs of SeqAIJ matrices

C = A * B

#include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/

#include <../src/mat/utils/freespace.h>

#include <petscbt.h>

#include <petsc/private/isimpl.h>

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

PETSC_INTERNextern __attribute__((visibility ("hidden"))) PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)

{

PetscErrorCode ierr;

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 18; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

if (scall == MAT_INITIAL_MATRIX) {

ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultSymbolic].active) ? (*PetscLogPLB)((MAT_MatMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),20,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),21,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultSymbolic].active) ? (*PetscLogPLE)((MAT_MatMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),22,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

}

ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultNumeric].active) ? (*PetscLogPLB)((MAT_MatMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),25,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),26,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultNumeric].active) ? (*PetscLogPLE)((MAT_MatMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),27,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

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);

}

PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)

{

PetscErrorCode ierr;

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 35; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

if (C->ops->matmultnumeric) {

ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),37,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

} else {

ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(A,B,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),39,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

}

PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)

{

PetscErrorCode ierr;

#if !defined(PETSC_HAVE_HYPRE)

const char *algTypes[8] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge"};

PetscInt nalg = 8;

#else

const char *algTypes[9] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge","hypre"};

PetscInt nalg = 9;

#endif

PetscInt alg = 0; /* set default algorithm */

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 56; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatMatMult","Mat")0; do { PetscOptionItems PetscOptionsObjectBase; PetscOptionItems
*PetscOptionsObject = &PetscOptionsObjectBase; PetscMemzero
(PetscOptionsObject,sizeof(PetscOptionItems)); for (PetscOptionsObject
->count=(PetscOptionsPublish?-1:1); PetscOptionsObject->
count<2; PetscOptionsObject->count++) { PetscErrorCode _5_ierr
= PetscOptionsBegin_Private(PetscOptionsObject,PetscObjectComm
((PetscObject)A),((PetscObject)A)->prefix,"MatMatMult","Mat"
);do {if (__builtin_expect(!!(_5_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),57,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),57,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

ierr = PetscOptionsEList("-matmatmult_via","Algorithmic approach","MatMatMult",algTypes,nalg,algTypes[0],&alg,NULL)PetscOptionsEList_Private(PetscOptionsObject,"-matmatmult_via"
,"Algorithmic approach","MatMatMult",algTypes,nalg,algTypes[0
],&alg,((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),58,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

ierr = PetscOptionsEnd()_5_ierr = PetscOptionsEnd_Private(PetscOptionsObject);do {if (
__builtin_expect(!!(_5_ierr),0)) return PetscError(((MPI_Comm
)0x44000001),59,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);}} while (0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),59,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

switch (alg) {

case 1:

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),62,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

break;

case 2:

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),65,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

break;

case 3:

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),68,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

break;

case 4:

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),71,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

break;

case 5:

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),74,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

break;

case 6:

ierr = MatMatMult_SeqAIJ_SeqAIJ_Combined(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),77,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

break;

case 7:

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),80,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

break;

#if defined(PETSC_HAVE_HYPRE)

case 8:

ierr = MatMatMultSymbolic_AIJ_AIJ_wHYPRE(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),84,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

break;

#endif

default:

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),88,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

break;

}

PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(Mat A,Mat B,PetscReal fill,Mat *C)

{

PetscErrorCode ierr;

Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;

PetscInt *ai=a->i,*bi=b->i,*ci,*cj;

PetscInt am =A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;

100

PetscReal afill;

101

PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;

102

PetscTable ta;

103

PetscBT lnkbt;

104

PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);

105

106

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 106; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

107

/* Get ci and cj */

108

/*---------------*/

109

/* Allocate ci array, arrays for fill computation and */

110

/* free space for accumulating nonzero column info */

111

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,111,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),111,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

112

ci[0] = 0;

113

114

/* create and initialize a linked list */

115

ierr = PetscTableCreate(bn,bn,&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),115,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

116

MatRowMergeMax_SeqAIJ(b,bm,ta){ PetscInt _j,_row,_nz,*_col; if (b) { for (_row=0; _row<bm
; _row++) { _nz = b->i[_row+1] - b->i[_row]; for (_j=0;
_j<_nz; _j++) { _col = _j + b->j + b->i[_row]; PetscTableAdd
(ta,*_col+1,1,INSERT_VALUES); } } } };

117

ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),117,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

118

ierr = PetscTableDestroy(&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),118,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

119

120

ierr = PetscLLCondensedCreate(Crmax,bn,&lnk,&lnkbt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),120,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

121

122

/* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */

123

ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),123,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

124

125

current_space = free_space;

126

127

/* Determine ci and cj */

128

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

129

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

130

aj = a->j + ai[i];

131

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

132

brow = aj[j];

133

bnzj = bi[brow+1] - bi[brow];

134

bj = b->j + bi[brow];

135

/* add non-zero cols of B into the sorted linked list lnk */

136

ierr = PetscLLCondensedAddSorted(bnzj,bj,lnk,lnkbt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),136,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

137

}

138

cnzi = lnk[0];

139

140

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

141

/* Double the amount of total space in the list */

142

if (current_space->local_remaining<cnzi) {

143

ierr = PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),143,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

144

ndouble++;

145

}

146

147

/* Copy data into free space, then initialize lnk */

148

ierr = PetscLLCondensedClean(bn,cnzi,current_space->array,lnk,lnkbt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),148,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

149

150

current_space->array += cnzi;

151

current_space->local_used += cnzi;

152

current_space->local_remaining -= cnzi;

153

154

ci[i+1] = ci[i] + cnzi;

155

}

156

157

/* Column indices are in the list of free space */

158

/* Allocate space for cj, initialize cj, and */

159

/* destroy list of free space and other temporary array(s) */

160

ierr = PetscMalloc1(ci[am]+1,&cj)PetscMallocA(1,PETSC_FALSE,160,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&cj)),(&cj));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),160,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

161

ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),161,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

162

ierr = PetscLLCondensedDestroy(lnk,lnkbt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),162,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

163

164

/* put together the new symbolic matrix */

165

ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),165,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

166

ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),166,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

167

ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),167,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

168

169

/* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */

170

/* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */

171

c = (Mat_SeqAIJ*)((*C)->data);

172

c->free_a = PETSC_FALSE;

173

c->free_ij = PETSC_TRUE;

174

c->nonew = 0;

175

(*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted; /* fast, needs non-scalable O(bn) array 'abdense' */

176

177

/* set MatInfo */

178

afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;

179

if (afill < 1.0) afill = 1.0;

180

c->maxnz = ci[am];

181

c->nz = ci[am];

182

(*C)->info.mallocs = ndouble;

183

(*C)->info.fill_ratio_given = fill;

184

(*C)->info.fill_ratio_needed = afill;

185

186

#if defined(PETSC_USE_INFO1)

187

if (ci[am]) {

188

ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),188,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

189

ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),189,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

190

} else {

191

ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),191,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

192

}

193

#endif

194

195

}

196

197

PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,Mat C)

198

{

199

PetscErrorCode ierr;

200

PetscLogDouble flops=0.0;

201

Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;

202

Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;

203

Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;

204

PetscInt *ai =a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;

205

PetscInt am =A->rmap->n,cm=C->rmap->n;

206

PetscInt i,j,k,anzi,bnzi,cnzi,brow;

207

PetscScalar *aa=a->a,*ba=b->a,*baj,*ca,valtmp;

208

PetscScalar *ab_dense;

209

210

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 210; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

211

if (!c->a) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */

212

ierr = PetscMalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_FALSE,212,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[cm]+1)*sizeof(**(&ca)),(&ca));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),212,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

213

c->a = ca;

214

c->free_a = PETSC_TRUE;

215

} else {

216

ca = c->a;

217

}

218

if (!c->matmult_abdense) {

219

ierr = PetscCalloc1(B->cmap->N,&ab_dense)PetscMallocA(1,PETSC_TRUE,219,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(B->cmap->N)*sizeof(**(&ab_dense)),(&ab_dense
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),219,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

220

c->matmult_abdense = ab_dense;

221

} else {

222

ab_dense = c->matmult_abdense;

223

}

224

225

/* clean old values in C */

226

ierr = PetscArrayzero(ca,ci[cm])PetscMemzero(ca,(ci[cm])*sizeof(*(ca)));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),226,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

227

/* Traverse A row-wise. */

228

/* Build the ith row in C by summing over nonzero columns in A, */

229

/* the rows of B corresponding to nonzeros of A. */

230

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

231

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

232

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

233

brow = aj[j];

234

bnzi = bi[brow+1] - bi[brow];

235

bjj = bj + bi[brow];

236

baj = ba + bi[brow];

237

/* perform dense axpy */

238

valtmp = aa[j];

239

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

240

ab_dense[bjj[k]] += valtmp*baj[k];

241

}

242

flops += 2*bnzi;

243

}

244

aj += anzi; aa += anzi;

245

246

cnzi = ci[i+1] - ci[i];

247

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

248

ca[k] += ab_dense[cj[k]];

249

ab_dense[cj[k]] = 0.0; /* zero ab_dense */

250

}

251

flops += cnzi;

252

cj += cnzi; ca += cnzi;

253

}

254

ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),254,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

255

ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),255,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

256

ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),256,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

257

258

}

259

260

PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,Mat C)

261

{

262

PetscErrorCode ierr;

263

PetscLogDouble flops=0.0;

264

Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;

265

Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;

266

Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;

267

PetscInt *ai = a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;

268

PetscInt am = A->rmap->N,cm=C->rmap->N;

269

PetscInt i,j,k,anzi,bnzi,cnzi,brow;

270

PetscScalar *aa=a->a,*ba=b->a,*baj,*ca=c->a,valtmp;

271

PetscInt nextb;

272

273

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 273; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

274

if (!ca) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */

275

ierr = PetscMalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_FALSE,275,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[cm]+1)*sizeof(**(&ca)),(&ca));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),275,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

276

c->a = ca;

277

c->free_a = PETSC_TRUE;

278

}

279

280

/* clean old values in C */

281

ierr = PetscArrayzero(ca,ci[cm])PetscMemzero(ca,(ci[cm])*sizeof(*(ca)));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),281,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

282

/* Traverse A row-wise. */

283

/* Build the ith row in C by summing over nonzero columns in A, */

284

/* the rows of B corresponding to nonzeros of A. */

285

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

286

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

287

cnzi = ci[i+1] - ci[i];

288

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

289

brow = aj[j];

290

bnzi = bi[brow+1] - bi[brow];

291

bjj = bj + bi[brow];

292

baj = ba + bi[brow];

293

/* perform sparse axpy */

294

valtmp = aa[j];

295

nextb = 0;

296

for (k=0; nextb<bnzi; k++) {

297

if (cj[k] == bjj[nextb]) { /* ccol == bcol */

298

ca[k] += valtmp*baj[nextb++];

299

}

300

}

301

flops += 2*bnzi;

302

}

303

aj += anzi; aa += anzi;

304

cj += cnzi; ca += cnzi;

305

}

306

307

ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),307,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

308

ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),308,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

309

ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),309,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

310

311

}

312

313

PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(Mat A,Mat B,PetscReal fill,Mat *C)

314

{

315

PetscErrorCode ierr;

316

Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;

317

PetscInt *ai = a->i,*bi=b->i,*ci,*cj;

318

PetscInt am = A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;

319

MatScalar *ca;

320

PetscReal afill;

321

PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;

322

PetscTable ta;

323

PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);

324

325

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 325; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

326

/* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_fast() */

327

/*-----------------------------------------------------------------------------------------*/

328

/* Allocate arrays for fill computation and free space for accumulating nonzero column */

329

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,329,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),329,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

330

ci[0] = 0;

331

332

/* create and initialize a linked list */

333

ierr = PetscTableCreate(bn,bn,&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),333,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

334

335

ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),335,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

336

ierr = PetscTableDestroy(&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),336,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

337

338

ierr = PetscLLCondensedCreate_fast(Crmax,&lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),338,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

339

340

/* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */

341

ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),341,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

342

current_space = free_space;

343

344

/* Determine ci and cj */

345

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

346

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

347

aj = a->j + ai[i];

348

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

349

brow = aj[j];

350

bnzj = bi[brow+1] - bi[brow];

351

bj = b->j + bi[brow];

352

/* add non-zero cols of B into the sorted linked list lnk */

353

ierr = PetscLLCondensedAddSorted_fast(bnzj,bj,lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),353,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

354

}

355

cnzi = lnk[1];

356

357

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

358

/* Double the amount of total space in the list */

359

if (current_space->local_remaining<cnzi) {

360

ierr = PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),360,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

361

ndouble++;

362

}

363

364

/* Copy data into free space, then initialize lnk */

365

ierr = PetscLLCondensedClean_fast(cnzi,current_space->array,lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),365,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

366

367

current_space->array += cnzi;

368

current_space->local_used += cnzi;

369

current_space->local_remaining -= cnzi;

370

371

ci[i+1] = ci[i] + cnzi;

372

}

373

374

/* Column indices are in the list of free space */

375

/* Allocate space for cj, initialize cj, and */

376

/* destroy list of free space and other temporary array(s) */

377

ierr = PetscMalloc1(ci[am]+1,&cj)PetscMallocA(1,PETSC_FALSE,377,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&cj)),(&cj));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),377,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

378

ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),378,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

379

ierr = PetscLLCondensedDestroy_fast(lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),379,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

380

381

/* Allocate space for ca */

382

ierr = PetscCalloc1(ci[am]+1,&ca)PetscMallocA(1,PETSC_TRUE,382,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&ca)),(&ca));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),382,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

383

384

/* put together the new symbolic matrix */

385

ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,ca,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),385,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

386

ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),386,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

387

ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),387,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

388

389

/* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */

390

/* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */

391

c = (Mat_SeqAIJ*)((*C)->data);

392

c->free_a = PETSC_TRUE;

393

c->free_ij = PETSC_TRUE;

394

c->nonew = 0;

395

396

(*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */

397

398

/* set MatInfo */

399

afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;

400

if (afill < 1.0) afill = 1.0;

401

c->maxnz = ci[am];

402

c->nz = ci[am];

403

(*C)->info.mallocs = ndouble;

404

(*C)->info.fill_ratio_given = fill;

405

(*C)->info.fill_ratio_needed = afill;

406

407

#if defined(PETSC_USE_INFO1)

408

if (ci[am]) {

409

ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),409,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

410

ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),410,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

411

} else {

412

ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),412,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

413

}

414

#endif

415

416

}

417

418

419

PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,PetscReal fill,Mat *C)

420

{

421

PetscErrorCode ierr;

422

Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;

423

PetscInt *ai = a->i,*bi=b->i,*ci,*cj;

424

PetscInt am = A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;

425

MatScalar *ca;

426

PetscReal afill;

427

PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;

428

PetscTable ta;

429

PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);

430

431

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 431; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

432

/* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_Scalalbe() */

433

/*---------------------------------------------------------------------------------------------*/

434

/* Allocate arrays for fill computation and free space for accumulating nonzero column */

435

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,435,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),435,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

436

ci[0] = 0;

437

438

/* create and initialize a linked list */

439

ierr = PetscTableCreate(bn,bn,&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),439,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

440

441

ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),441,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

442

ierr = PetscTableDestroy(&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),442,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

443

ierr = PetscLLCondensedCreate_Scalable(Crmax,&lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),443,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

444

445

/* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */

446

ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),446,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

447

current_space = free_space;

448

449

/* Determine ci and cj */

450

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

451

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

452

aj = a->j + ai[i];

453

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

454

brow = aj[j];

455

bnzj = bi[brow+1] - bi[brow];

456

bj = b->j + bi[brow];

457

/* add non-zero cols of B into the sorted linked list lnk */

458

ierr = PetscLLCondensedAddSorted_Scalable(bnzj,bj,lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),458,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

459

}

460

cnzi = lnk[0];

461

462

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

463

/* Double the amount of total space in the list */

464

if (current_space->local_remaining<cnzi) {

465

ierr = PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),465,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

466

ndouble++;

467

}

468

469

/* Copy data into free space, then initialize lnk */

470

ierr = PetscLLCondensedClean_Scalable(cnzi,current_space->array,lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),470,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

471

472

current_space->array += cnzi;

473

current_space->local_used += cnzi;

474

current_space->local_remaining -= cnzi;

475

476

ci[i+1] = ci[i] + cnzi;

477

}

478

479

/* Column indices are in the list of free space */

480

/* Allocate space for cj, initialize cj, and */

481

/* destroy list of free space and other temporary array(s) */

482

ierr = PetscMalloc1(ci[am]+1,&cj)PetscMallocA(1,PETSC_FALSE,482,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&cj)),(&cj));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),482,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

483

ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),483,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

484

ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),484,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

485

486

/* Allocate space for ca */

487

/*-----------------------*/

488

ierr = PetscCalloc1(ci[am]+1,&ca)PetscMallocA(1,PETSC_TRUE,488,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&ca)),(&ca));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),488,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

489

490

/* put together the new symbolic matrix */

491

ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,ca,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),491,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

492

ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),492,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

493

ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),493,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

494

495

/* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */

496

/* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */

497

c = (Mat_SeqAIJ*)((*C)->data);

498

c->free_a = PETSC_TRUE;

499

c->free_ij = PETSC_TRUE;

500

c->nonew = 0;

501

502

(*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */

503

504

/* set MatInfo */

505

afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;

506

if (afill < 1.0) afill = 1.0;

507

c->maxnz = ci[am];

508

c->nz = ci[am];

509

(*C)->info.mallocs = ndouble;

510

(*C)->info.fill_ratio_given = fill;

511

(*C)->info.fill_ratio_needed = afill;

512

513

#if defined(PETSC_USE_INFO1)

514

if (ci[am]) {

515

ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),515,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

516

ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),516,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

517

} else {

518

ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),518,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

519

}

520

#endif

521

522

}

523

524

PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(Mat A,Mat B,PetscReal fill,Mat *C)

525

{

526

PetscErrorCode ierr;

527

Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;

528

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

529

PetscInt *ci,*cj,*bb;

530

PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;

531

PetscReal afill;

532

PetscInt i,j,col,ndouble = 0;

533

PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);

534

PetscHeap h;

535

536

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 536; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

537

/* Get ci and cj - by merging sorted rows using a heap */

538

/*---------------------------------------------------------------------------------------------*/

539

/* Allocate arrays for fill computation and free space for accumulating nonzero column */

540

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,540,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),540,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

541

ci[0] = 0;

542

543

/* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */

544

ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),544,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

545

current_space = free_space;

546

547

ierr = PetscHeapCreate(a->rmax,&h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),547,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

548

ierr = PetscMalloc1(a->rmax,&bb)PetscMallocA(1,PETSC_FALSE,548,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a->rmax)*sizeof(**(&bb)),(&bb));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),548,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

549

550

/* Determine ci and cj */

551

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

552

const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */

553

const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */

554

ci[i+1] = ci[i];

555

/* Populate the min heap */

556

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

557

bb[j] = bi[acol[j]]; /* bb points at the start of the row */

558

if (bb[j] < bi[acol[j]+1]) { /* Add if row is nonempty */

559

ierr = PetscHeapAdd(h,j,bj[bb[j]++]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),559,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

560

}

561

}

562

/* Pick off the min element, adding it to free space */

563

ierr = PetscHeapPop(h,&j,&col);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),563,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

564

while (j >= 0) {

565

if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */

566

ierr = PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20)(((PetscIntMultTruncate(2,current_space->total_array_size)
)<(16 << 20)) ? (PetscIntMultTruncate(2,current_space
->total_array_size)) : (16 << 20)),&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),566,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

567

ndouble++;

568

}

569

*(current_space->array++) = col;

570

current_space->local_used++;

571

current_space->local_remaining--;

572

ci[i+1]++;

573

574

/* stash if anything else remains in this row of B */

575

if (bb[j] < bi[acol[j]+1]) {ierr = PetscHeapStash(h,j,bj[bb[j]++]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),575,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);}

576

while (1) { /* pop and stash any other rows of B that also had an entry in this column */

577

PetscInt j2,col2;

578

ierr = PetscHeapPeek(h,&j2,&col2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),578,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

579

if (col2 != col) break;

580

ierr = PetscHeapPop(h,&j2,&col2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),580,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

581

if (bb[j2] < bi[acol[j2]+1]) {ierr = PetscHeapStash(h,j2,bj[bb[j2]++]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),581,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);}

582

}

583

/* Put any stashed elements back into the min heap */

584

ierr = PetscHeapUnstash(h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),584,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

585

ierr = PetscHeapPop(h,&j,&col);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),585,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

586

}

587

}

588

ierr = PetscFree(bb)((*PetscTrFree)((void*)(bb),588,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((bb) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),588,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

589

ierr = PetscHeapDestroy(&h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),589,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

590

591

/* Column indices are in the list of free space */

592

/* Allocate space for cj, initialize cj, and */

593

/* destroy list of free space and other temporary array(s) */

594

ierr = PetscMalloc1(ci[am],&cj)PetscMallocA(1,PETSC_FALSE,594,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am])*sizeof(**(&cj)),(&cj));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),594,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

595

ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),595,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

596

597

/* put together the new symbolic matrix */

598

ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),598,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

599

ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),599,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

600

ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),600,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

601

602

/* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */

603

/* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */

604

c = (Mat_SeqAIJ*)((*C)->data);

605

c->free_a = PETSC_TRUE;

606

c->free_ij = PETSC_TRUE;

607

c->nonew = 0;

608

609

(*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

610

611

/* set MatInfo */

612

afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;

613

if (afill < 1.0) afill = 1.0;

614

c->maxnz = ci[am];

615

c->nz = ci[am];

616

(*C)->info.mallocs = ndouble;

617

(*C)->info.fill_ratio_given = fill;

618

(*C)->info.fill_ratio_needed = afill;

619

620

#if defined(PETSC_USE_INFO1)

621

if (ci[am]) {

622

ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),622,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

623

ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),623,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

624

} else {

625

ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),625,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

626

}

627

#endif

628

629

}

630

631

PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(Mat A,Mat B,PetscReal fill,Mat *C)

632

{

633

PetscErrorCode ierr;

634

Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;

635

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

636

PetscInt *ci,*cj,*bb;

637

PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;

638

PetscReal afill;

639

PetscInt i,j,col,ndouble = 0;

640

PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);

641

PetscHeap h;

642

PetscBT bt;

643

644

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 644; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

645

/* Get ci and cj - using a heap for the sorted rows, but use BT so that each index is only added once */

646

/*---------------------------------------------------------------------------------------------*/

647

/* Allocate arrays for fill computation and free space for accumulating nonzero column */

648

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,648,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),648,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

649

ci[0] = 0;

650

651

/* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */

652

ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),652,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

653

654

current_space = free_space;

Value assigned to 'current_space'

→

655

656

ierr = PetscHeapCreate(a->rmax,&h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),656,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

657

ierr = PetscMalloc1(a->rmax,&bb)PetscMallocA(1,PETSC_FALSE,657,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a->rmax)*sizeof(**(&bb)),(&bb));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),657,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

658

ierr = PetscBTCreate(bn,&bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),658,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

659

660

/* Determine ci and cj */

661

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

←

Assuming 'i' is < 'am'

→

←

Loop condition is true. Entering loop body

→

←

Assuming 'i' is < 'am'

→

←

Loop condition is true. Entering loop body

→

662

const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */

663

const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */

664

const PetscInt *fptr = current_space->array; /* Save beginning of the row so we can clear the BT later */

←

Assuming pointer value is null

→

665

ci[i+1] = ci[i];

666

/* Populate the min heap */

667

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

←

Assuming 'j' is >= 'anzi'

→

←

Loop condition is false. Execution continues on line 679

→

←

Assuming 'j' is >= 'anzi'

→

←

Loop condition is false. Execution continues on line 679

→

668

PetscInt brow = acol[j];

669

for (bb[j] = bi[brow]; bb[j] < bi[brow+1]; bb[j]++) {

670

PetscInt bcol = bj[bb[j]];

671

if (!PetscBTLookupSet(bt,bcol)) { /* new entry */

672

ierr = PetscHeapAdd(h,j,bcol);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),672,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

673

bb[j]++;

674

break;

675

}

676

}

677

}

678

/* Pick off the min element, adding it to free space */

679

ierr = PetscHeapPop(h,&j,&col);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),679,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

680

while (j >= 0) {

←

Assuming 'j' is < 0

→

←

Loop condition is false. Execution continues on line 702

Assuming 'j' is >= 0

Loop condition is true. Entering loop body

→

681

if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */

←

Assuming the condition is false

→

←

Taking false branch

→

682

fptr = NULL((void*)0); /* need PetscBTMemzero */

683

ierr = PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20)(((PetscIntMultTruncate(2,current_space->total_array_size)
)<(16 << 20)) ? (PetscIntMultTruncate(2,current_space
->total_array_size)) : (16 << 20)),&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),683,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

684

ndouble++;

685

}

686

*(current_space->array++) = col;

←

Dereference of null pointer

687

current_space->local_used++;

688

current_space->local_remaining--;

689

ci[i+1]++;

690

691

/* stash if anything else remains in this row of B */

692

for (; bb[j] < bi[acol[j]+1]; bb[j]++) {

693

PetscInt bcol = bj[bb[j]];

694

if (!PetscBTLookupSet(bt,bcol)) { /* new entry */

695

ierr = PetscHeapAdd(h,j,bcol);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),695,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

696

bb[j]++;

697

break;

698

}

699

}

700

ierr = PetscHeapPop(h,&j,&col);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),700,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

701

}

702

if (fptr) { /* Clear the bits for this row */

←

Assuming 'fptr' is null

→

←

Taking false branch

→

703

for (; fptr<current_space->array; fptr++) {ierr = PetscBTClear(bt,*fptr);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),703,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);}

704

} else { /* We reallocated so we don't remember (easily) how to clear only the bits we changed */

705

ierr = PetscBTMemzero(bn,bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),705,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

706

}

707

}

708

ierr = PetscFree(bb)((*PetscTrFree)((void*)(bb),708,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((bb) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),708,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

709

ierr = PetscHeapDestroy(&h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),709,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

710

ierr = PetscBTDestroy(&bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),710,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

711

712

/* Column indices are in the list of free space */

713

/* Allocate space for cj, initialize cj, and */

714

/* destroy list of free space and other temporary array(s) */

715

ierr = PetscMalloc1(ci[am],&cj)PetscMallocA(1,PETSC_FALSE,715,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am])*sizeof(**(&cj)),(&cj));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),715,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

716

ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),716,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

717

718

/* put together the new symbolic matrix */

719

ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),719,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

720

ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),720,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

721

ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),721,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

722

723

/* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */

724

/* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */

725

c = (Mat_SeqAIJ*)((*C)->data);

726

c->free_a = PETSC_TRUE;

727

c->free_ij = PETSC_TRUE;

728

c->nonew = 0;

729

730

(*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

731

732

/* set MatInfo */

733

afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;

734

if (afill < 1.0) afill = 1.0;

735

c->maxnz = ci[am];

736

c->nz = ci[am];

737

(*C)->info.mallocs = ndouble;

738

(*C)->info.fill_ratio_given = fill;

739

(*C)->info.fill_ratio_needed = afill;

740

741

#if defined(PETSC_USE_INFO1)

742

if (ci[am]) {

743

ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),743,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

744

ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),744,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

745

} else {

746

ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),746,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

747

}

748

#endif

749

750

}

751

752

753

PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(Mat A,Mat B,PetscReal fill,Mat *C)

754

{

755

PetscErrorCode ierr;

756

Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;

757

const PetscInt *ai=a->i,*bi=b->i,*aj=a->j,*bj=b->j,*inputi,*inputj,*inputcol,*inputcol_L1;

758

PetscInt *ci,*cj,*outputj,worki_L1[9],worki_L2[9];

759

PetscInt c_maxmem,a_maxrownnz=0,a_rownnz;

760

const PetscInt workcol[8]={0,1,2,3,4,5,6,7};

761

const PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;

762

const PetscInt *brow_ptr[8],*brow_end[8];

763

PetscInt window[8];

764

PetscInt window_min,old_window_min,ci_nnz,outputi_nnz=0,L1_nrows,L2_nrows;

765

PetscInt i,k,ndouble=0,L1_rowsleft,rowsleft;

766

PetscReal afill;

767

PetscInt *workj_L1,*workj_L2,*workj_L3;

768

PetscInt L1_nnz,L2_nnz;

769

770

/* Step 1: Get upper bound on memory required for allocation.

771

Because of the way virtual memory works,

772

only the memory pages that are actually needed will be physically allocated. */

773

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 773; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

774

ierr = PetscMalloc1(am+1,&ci)PetscMallocA(1,PETSC_FALSE,774,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+1)*sizeof(**(&ci)),(&ci));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),774,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

775

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

776

const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */

777

const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */

778

a_rownnz = 0;

779

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

780

a_rownnz += bi[acol[k]+1] - bi[acol[k]];

781

if (a_rownnz > bn) {

782

a_rownnz = bn;

783

break;

784

}

785

}

786

a_maxrownnz = PetscMax(a_maxrownnz, a_rownnz)(((a_maxrownnz)<(a_rownnz)) ? (a_rownnz) : (a_maxrownnz));

787

}

788

/* temporary work areas for merging rows */

789

ierr = PetscMalloc1(a_maxrownnz*8,&workj_L1)PetscMallocA(1,PETSC_FALSE,789,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a_maxrownnz*8)*sizeof(**(&workj_L1)),(&workj_L1
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),789,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

790

ierr = PetscMalloc1(a_maxrownnz*8,&workj_L2)PetscMallocA(1,PETSC_FALSE,790,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a_maxrownnz*8)*sizeof(**(&workj_L2)),(&workj_L2
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),790,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

791

ierr = PetscMalloc1(a_maxrownnz,&workj_L3)PetscMallocA(1,PETSC_FALSE,791,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a_maxrownnz)*sizeof(**(&workj_L3)),(&workj_L3
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),791,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

792

793

/* This should be enough for almost all matrices. If not, memory is reallocated later. */

794

c_maxmem = 8*(ai[am]+bi[bm]);

795

/* Step 2: Populate pattern for C */

796

ierr = PetscMalloc1(c_maxmem,&cj)PetscMallocA(1,PETSC_FALSE,796,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(c_maxmem)*sizeof(**(&cj)),(&cj));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),796,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

797

798

ci_nnz = 0;

799

ci[0] = 0;

800

worki_L1[0] = 0;

801

worki_L2[0] = 0;

802

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

803

const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */

804

const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */

805

rowsleft = anzi;

806

inputcol_L1 = acol;

807

L2_nnz = 0;

808

L2_nrows = 1; /* Number of rows to be merged on Level 3. output of L3 already exists -> initial value 1 */

809

worki_L2[1] = 0;

810

outputi_nnz = 0;

811

812

/* If the number of indices in C so far + the max number of columns in the next row > c_maxmem -> allocate more memory */

813

while (ci_nnz+a_maxrownnz > c_maxmem) {

814

c_maxmem *= 2;

815

ndouble++;

816

ierr = PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj)((*PetscTrRealloc)((sizeof(PetscInt)*c_maxmem),816,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(void**)(&cj)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),816,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

817

}

818

819

while (rowsleft) {

820

L1_rowsleft = PetscMin(64, rowsleft)(((64)<(rowsleft)) ? (64) : (rowsleft)); /* In the inner loop max 64 rows of B can be merged */

821

L1_nrows = 0;

822

L1_nnz = 0;

823

inputcol = inputcol_L1;

824

inputi = bi;

825

inputj = bj;

826

827

/* The following macro is used to specialize for small rows in A.

828

This helps with compiler unrolling, improving performance substantially.

829

Input: inputj inputi inputcol bn

830

Output: outputj outputi_nnz */

831

#define MatMatMultSymbolic_RowMergeMacro(ANNZ) \

832

window_min = bn; \

833

outputi_nnz = 0; \

834

for (k=0; k<ANNZ; ++k) { \

835

brow_ptr[k] = inputj + inputi[inputcol[k]]; \

836

brow_end[k] = inputj + inputi[inputcol[k]+1]; \

837

window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \

838

window_min = PetscMin(window[k], window_min)(((window[k])<(window_min)) ? (window[k]) : (window_min)); \

839

} \

840

while (window_min < bn) { \

841

outputj[outputi_nnz++] = window_min; \

842

/* advance front and compute new minimum */ \

843

old_window_min = window_min; \

844

window_min = bn; \

845

for (k=0; k<ANNZ; ++k) { \

846

if (window[k] == old_window_min) { \

847

brow_ptr[k]++; \

848

window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \

849

} \

850

window_min = PetscMin(window[k], window_min)(((window[k])<(window_min)) ? (window[k]) : (window_min)); \

851

} \

852

}

853

854

/************** L E V E L 1 ***************/

855

/* Merge up to 8 rows of B to L1 work array*/

856

while (L1_rowsleft) {

857

outputi_nnz = 0;

858

if (anzi > 8) outputj = workj_L1 + L1_nnz; /* Level 1 rowmerge*/

859

else outputj = cj + ci_nnz; /* Merge directly to C */

860

861

switch (L1_rowsleft) {

862

case 1: brow_ptr[0] = inputj + inputi[inputcol[0]];

863

brow_end[0] = inputj + inputi[inputcol[0]+1];

864

for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */

865

inputcol += L1_rowsleft;

866

rowsleft -= L1_rowsleft;

867

L1_rowsleft = 0;

868

break;

869

case 2: MatMatMultSymbolic_RowMergeMacro(2);

870

inputcol += L1_rowsleft;

871

rowsleft -= L1_rowsleft;

872

L1_rowsleft = 0;

873

break;

874

case 3: MatMatMultSymbolic_RowMergeMacro(3);

875

inputcol += L1_rowsleft;

876

rowsleft -= L1_rowsleft;

877

L1_rowsleft = 0;

878

break;

879

case 4: MatMatMultSymbolic_RowMergeMacro(4);

880

inputcol += L1_rowsleft;

881

rowsleft -= L1_rowsleft;

882

L1_rowsleft = 0;

883

break;

884

case 5: MatMatMultSymbolic_RowMergeMacro(5);

885

inputcol += L1_rowsleft;

886

rowsleft -= L1_rowsleft;

887

L1_rowsleft = 0;

888

break;

889

case 6: MatMatMultSymbolic_RowMergeMacro(6);

890

inputcol += L1_rowsleft;

891

rowsleft -= L1_rowsleft;

892

L1_rowsleft = 0;

893

break;

894

case 7: MatMatMultSymbolic_RowMergeMacro(7);

895

inputcol += L1_rowsleft;

896

rowsleft -= L1_rowsleft;

897

L1_rowsleft = 0;

898

break;

899

default: MatMatMultSymbolic_RowMergeMacro(8);

900

inputcol += 8;

901

rowsleft -= 8;

902

L1_rowsleft -= 8;

903

break;

904

}

905

inputcol_L1 = inputcol;

906

L1_nnz += outputi_nnz;

907

worki_L1[++L1_nrows] = L1_nnz;

908

}

909

910

/********************** L E V E L 2 ************************/

911

/* Merge from L1 work array to either C or to L2 work array */

912

if (anzi > 8) {

913

inputi = worki_L1;

914

inputj = workj_L1;

915

inputcol = workcol;

916

outputi_nnz = 0;

917

918

if (anzi <= 64) outputj = cj + ci_nnz; /* Merge from L1 work array to C */

919

else outputj = workj_L2 + L2_nnz; /* Merge from L1 work array to L2 work array */

920

921

switch (L1_nrows) {

922

case 1: brow_ptr[0] = inputj + inputi[inputcol[0]];

923

brow_end[0] = inputj + inputi[inputcol[0]+1];

924

for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */

925

break;

926

case 2: MatMatMultSymbolic_RowMergeMacro(2); break;

927

case 3: MatMatMultSymbolic_RowMergeMacro(3); break;

928

case 4: MatMatMultSymbolic_RowMergeMacro(4); break;

929

case 5: MatMatMultSymbolic_RowMergeMacro(5); break;

930

case 6: MatMatMultSymbolic_RowMergeMacro(6); break;

931

case 7: MatMatMultSymbolic_RowMergeMacro(7); break;

932

case 8: MatMatMultSymbolic_RowMergeMacro(8); break;

933

default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatMult logic error: Not merging 1-8 rows from L1 work array!")return PetscError(((MPI_Comm)0x44000001),933,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,56,PETSC_ERROR_INITIAL,"MatMatMult logic error: Not merging 1-8 rows from L1 work array!"
);

934

}

935

L2_nnz += outputi_nnz;

936

worki_L2[++L2_nrows] = L2_nnz;

937

938

/************************ L E V E L 3 **********************/

939

/* Merge from L2 work array to either C or to L2 work array */

940

if (anzi > 64 && (L2_nrows == 8 || rowsleft == 0)) {

941

inputi = worki_L2;

942

inputj = workj_L2;

943

inputcol = workcol;

944

outputi_nnz = 0;

945

if (rowsleft) outputj = workj_L3;

946

else outputj = cj + ci_nnz;

947

switch (L2_nrows) {

948

case 1: brow_ptr[0] = inputj + inputi[inputcol[0]];

949

brow_end[0] = inputj + inputi[inputcol[0]+1];

950

for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */

951

break;

952

case 2: MatMatMultSymbolic_RowMergeMacro(2); break;

953

case 3: MatMatMultSymbolic_RowMergeMacro(3); break;

954

case 4: MatMatMultSymbolic_RowMergeMacro(4); break;

955

case 5: MatMatMultSymbolic_RowMergeMacro(5); break;

956

case 6: MatMatMultSymbolic_RowMergeMacro(6); break;

957

case 7: MatMatMultSymbolic_RowMergeMacro(7); break;

958

case 8: MatMatMultSymbolic_RowMergeMacro(8); break;

959

default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatMult logic error: Not merging 1-8 rows from L2 work array!")return PetscError(((MPI_Comm)0x44000001),959,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,56,PETSC_ERROR_INITIAL,"MatMatMult logic error: Not merging 1-8 rows from L2 work array!"
);

960

}

961

L2_nrows = 1;

962

L2_nnz = outputi_nnz;

963

worki_L2[1] = outputi_nnz;

964

/* Copy to workj_L2 */

965

if (rowsleft) {

966

for (k=0; k<outputi_nnz; ++k) workj_L2[k] = outputj[k];

967

}

968

}

969

}

970

} /* while (rowsleft) */

971

#undef MatMatMultSymbolic_RowMergeMacro

972

973

/* terminate current row */

974

ci_nnz += outputi_nnz;

975

ci[i+1] = ci_nnz;

976

}

977

978

/* Step 3: Create the new symbolic matrix */

979

ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),979,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

980

ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),980,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

981

ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),981,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

982

983

/* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */

984

/* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */

985

c = (Mat_SeqAIJ*)((*C)->data);

986

c->free_a = PETSC_TRUE;

987

c->free_ij = PETSC_TRUE;

988

c->nonew = 0;

989

990

(*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

991

992

/* set MatInfo */

993

afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;

994

if (afill < 1.0) afill = 1.0;

995

c->maxnz = ci[am];

996

c->nz = ci[am];

997

(*C)->info.mallocs = ndouble;

998

(*C)->info.fill_ratio_given = fill;

999

(*C)->info.fill_ratio_needed = afill;

1000

1001

#if defined(PETSC_USE_INFO1)

1002

if (ci[am]) {

1003

ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1003,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1004

ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1004,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1005

} else {

1006

ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1006,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1007

}

1008

#endif

1009

1010

/* Step 4: Free temporary work areas */

1011

ierr = PetscFree(workj_L1)((*PetscTrFree)((void*)(workj_L1),1011,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((workj_L1) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1011,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1012

ierr = PetscFree(workj_L2)((*PetscTrFree)((void*)(workj_L2),1012,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((workj_L2) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1012,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1013

ierr = PetscFree(workj_L3)((*PetscTrFree)((void*)(workj_L3),1013,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((workj_L3) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1013,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1014

1015

}

1016

1017

/* concatenate unique entries and then sort */

1018

PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,PetscReal fill,Mat *C)

1019

{

1020

PetscErrorCode ierr;

1021

Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;

1022

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

1023

PetscInt *ci,*cj;

1024

PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;

1025

PetscReal afill;

1026

PetscInt i,j,ndouble = 0;

1027

PetscSegBuffer seg,segrow;

1028

char *seen;

1029

1030

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1030; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1031

ierr = PetscMalloc1(am+1,&ci)PetscMallocA(1,PETSC_FALSE,1031,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+1)*sizeof(**(&ci)),(&ci));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1031,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1032

ci[0] = 0;

1033

1034

/* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */

1035

ierr = PetscSegBufferCreate(sizeof(PetscInt),(PetscInt)(fill*(ai[am]+bi[bm])),&seg);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1035,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1036

ierr = PetscSegBufferCreate(sizeof(PetscInt),100,&segrow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1036,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1037

ierr = PetscCalloc1(bn,&seen)PetscMallocA(1,PETSC_TRUE,1037,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(bn)*sizeof(**(&seen)),(&seen));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1037,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1038

1039

/* Determine ci and cj */

1040

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

1041

const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */

1042

const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */

1043

PetscInt packlen = 0,*PETSC_RESTRICT__restrict crow;

1044

/* Pack segrow */

1045

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

1046

PetscInt brow = acol[j],bjstart = bi[brow],bjend = bi[brow+1],k;

1047

for (k=bjstart; k<bjend; k++) {

1048

PetscInt bcol = bj[k];

1049

if (!seen[bcol]) { /* new entry */

1050

PetscInt *PETSC_RESTRICT__restrict slot;

1051

ierr = PetscSegBufferGetInts(segrow,1,&slot);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1051,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1052

*slot = bcol;

1053

seen[bcol] = 1;

1054

packlen++;

1055

}

1056

}

1057

}

1058

ierr = PetscSegBufferGetInts(seg,packlen,&crow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1058,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1059

ierr = PetscSegBufferExtractTo(segrow,crow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1059,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1060

ierr = PetscSortInt(packlen,crow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1060,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1061

ci[i+1] = ci[i] + packlen;

1062

for (j=0; j<packlen; j++) seen[crow[j]] = 0;

1063

}

1064

ierr = PetscSegBufferDestroy(&segrow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1064,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1065

ierr = PetscFree(seen)((*PetscTrFree)((void*)(seen),1065,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((seen) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1065,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1066

1067

/* Column indices are in the segmented buffer */

1068

ierr = PetscSegBufferExtractAlloc(seg,&cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1068,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1069

ierr = PetscSegBufferDestroy(&seg);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1069,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1070

1071

/* put together the new symbolic matrix */

1072

ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1072,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1073

ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1073,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1074

ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1074,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1075

1076

/* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */

1077

/* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */

1078

c = (Mat_SeqAIJ*)((*C)->data);

1079

c->free_a = PETSC_TRUE;

1080

c->free_ij = PETSC_TRUE;

1081

c->nonew = 0;

1082

1083

(*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;

1084

1085

/* set MatInfo */

1086

afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;

1087

if (afill < 1.0) afill = 1.0;

1088

c->maxnz = ci[am];

1089

c->nz = ci[am];

1090

(*C)->info.mallocs = ndouble;

1091

(*C)->info.fill_ratio_given = fill;

1092

(*C)->info.fill_ratio_needed = afill;

1093

1094

#if defined(PETSC_USE_INFO1)

1095

if (ci[am]) {

1096

ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1096,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1097

ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1097,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1098

} else {

1099

ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1099,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1100

}

1101

#endif

1102

1103

}

1104

1105

/* This routine is not used. Should be removed! */

1106

PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)

1107

{

1108

PetscErrorCode ierr;

1109

1110

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1110; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1111

if (scall == MAT_INITIAL_MATRIX) {

1112

ierr = PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatTransposeMultSymbolic].active) ? (*PetscLogPLB)((MAT_MatTransposeMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1112,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1113

ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1113,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1114

ierr = PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatTransposeMultSymbolic].active) ? (*PetscLogPLE)((MAT_MatTransposeMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1114,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1115

}

1116

ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatTransposeMultNumeric].active) ? (*PetscLogPLB)((MAT_MatTransposeMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1116,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1117

ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1117,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1118

ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatTransposeMultNumeric].active) ? (*PetscLogPLE)((MAT_MatTransposeMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1118,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1119

1120

}

1121

1122

PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(Mat A)

1123

{

1124

PetscErrorCode ierr;

1125

Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;

1126

Mat_MatMatTransMult *abt=a->abt;

1127

1128

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1128; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1129

ierr = (abt->destroy)(A);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1129,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1130

ierr = MatTransposeColoringDestroy(&abt->matcoloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1130,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1131

ierr = MatDestroy(&abt->Bt_den);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1131,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1132

ierr = MatDestroy(&abt->ABt_den);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1132,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1133

ierr = PetscFree(abt)((*PetscTrFree)((void*)(abt),1133,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((abt) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1133,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1134

1135

}

1136

1137

PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)

1138

{

1139

PetscErrorCode ierr;

1140

Mat Bt;

1141

PetscInt *bti,*btj;

1142

Mat_MatMatTransMult *abt;

1143

Mat_SeqAIJ *c;

1144

1145

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1145; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1146

/* create symbolic Bt */

1147

ierr = MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1147,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1148

ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF((MPI_Comm)0x44000001),B->cmap->n,B->rmap->n,bti,btj,NULL((void*)0),&Bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1148,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1149

ierr = MatSetBlockSizes(Bt,PetscAbs(A->cmap->bs)(((A->cmap->bs) >= 0) ? (A->cmap->bs) : (-(A->
cmap->bs))),PetscAbs(B->cmap->bs)(((B->cmap->bs) >= 0) ? (B->cmap->bs) : (-(B->
cmap->bs))));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1149,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1150

ierr = MatSetType(Bt,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1150,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1151

1152

/* get symbolic C=A*Bt */

1153

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1153,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1154

1155

/* create a supporting struct for reuse intermidiate dense matrices with matcoloring */

1156

ierr = PetscNew(&abt)PetscMallocA(1,PETSC_TRUE,1156,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(1)*sizeof(**((&abt))),((&abt)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1156,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1157

c = (Mat_SeqAIJ*)(*C)->data;

1158

c->abt = abt;

1159

1160

abt->usecoloring = PETSC_FALSE;

1161

abt->destroy = (*C)->ops->destroy;

1162

(*C)->ops->destroy = MatDestroy_SeqAIJ_MatMatMultTrans;

1163

1164

ierr = PetscOptionsGetBool(((PetscObject)A)->options,NULL((void*)0),"-matmattransmult_color",&abt->usecoloring,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1164,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1165

if (abt->usecoloring) {

1166

/* Create MatTransposeColoring from symbolic C=A*B^T */

1167

MatTransposeColoring matcoloring;

1168

MatColoring coloring;

1169

ISColoring iscoloring;

1170

Mat Bt_dense,C_dense;

1171

Mat_SeqAIJ *c=(Mat_SeqAIJ*)(*C)->data;

1172

/* inode causes memory problem, don't know why */

1173

if (c->inode.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MAT_USE_INODES is not supported. Use '-mat_no_inode'")return PetscError(((MPI_Comm)0x44000001),1173,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,56,PETSC_ERROR_INITIAL,"MAT_USE_INODES is not supported. Use '-mat_no_inode'"
);

1174

1175

ierr = MatColoringCreate(*C,&coloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1175,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1176

ierr = MatColoringSetDistance(coloring,2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1176,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1177

ierr = MatColoringSetType(coloring,MATCOLORINGSL"sl");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1177,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1178

ierr = MatColoringSetFromOptions(coloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1178,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1179

ierr = MatColoringApply(coloring,&iscoloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1179,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1180

ierr = MatColoringDestroy(&coloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1180,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1181

ierr = MatTransposeColoringCreate(*C,iscoloring,&matcoloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1181,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1182

1183

abt->matcoloring = matcoloring;

1184

1185

ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1185,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1186

1187

/* Create Bt_dense and C_dense = A*Bt_dense */

1188

ierr = MatCreate(PETSC_COMM_SELF((MPI_Comm)0x44000001),&Bt_dense);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1188,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1189

ierr = MatSetSizes(Bt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1189,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1190

ierr = MatSetType(Bt_dense,MATSEQDENSE"seqdense");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1190,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1191

ierr = MatSeqDenseSetPreallocation(Bt_dense,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1191,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1192

1193

Bt_dense->assembled = PETSC_TRUE;

1194

abt->Bt_den = Bt_dense;

1195

1196

ierr = MatCreate(PETSC_COMM_SELF((MPI_Comm)0x44000001),&C_dense);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1196,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1197

ierr = MatSetSizes(C_dense,A->rmap->n,matcoloring->ncolors,A->rmap->n,matcoloring->ncolors);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1197,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1198

ierr = MatSetType(C_dense,MATSEQDENSE"seqdense");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1198,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1199

ierr = MatSeqDenseSetPreallocation(C_dense,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1199,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1200

1201

Bt_dense->assembled = PETSC_TRUE;

1202

abt->ABt_den = C_dense;

1203

1204

#if defined(PETSC_USE_INFO1)

1205

{

1206

Mat_SeqAIJ *c = (Mat_SeqAIJ*)(*C)->data;

1207

ierr = PetscInfo7(*C,"Use coloring of C=A*B^T; B^T: %D %D, Bt_dense: %D,%D; Cnz %D / (cm*ncolors %D) = %g\n",B->cmap->n,B->rmap->n,Bt_dense->rmap->n,Bt_dense->cmap->n,c->nz,A->rmap->n*matcoloring->ncolors,(PetscReal)(c->nz)/(A->rmap->n*matcoloring->ncolors))PetscInfo_Private(__func__,*C,"Use coloring of C=A*B^T; B^T: %D %D, Bt_dense: %D,%D; Cnz %D / (cm*ncolors %D) = %g\n"
,B->cmap->n,B->rmap->n,Bt_dense->rmap->n,Bt_dense
->cmap->n,c->nz,A->rmap->n*matcoloring->ncolors
,(PetscReal)(c->nz)/(A->rmap->n*matcoloring->ncolors
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1207,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1208

}

1209

#endif

1210

}

1211

/* clean up */

1212

ierr = MatDestroy(&Bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1212,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1213

ierr = MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1213,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1214

1215

}

1216

1217

PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)

1218

{

1219

PetscErrorCode ierr;

1220

Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;

1221

PetscInt *ai =a->i,*aj=a->j,*bi=b->i,*bj=b->j,anzi,bnzj,nexta,nextb,*acol,*bcol,brow;

1222

PetscInt cm =C->rmap->n,*ci=c->i,*cj=c->j,i,j,cnzi,*ccol;

1223

PetscLogDouble flops=0.0;

1224

MatScalar *aa =a->a,*aval,*ba=b->a,*bval,*ca,*cval;

1225

Mat_MatMatTransMult *abt = c->abt;

1226

1227

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1227; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1228

/* clear old values in C */

1229

if (!c->a) {

1230

ierr = PetscCalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_TRUE,1230,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[cm]+1)*sizeof(**(&ca)),(&ca));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1230,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1231

c->a = ca;

1232

c->free_a = PETSC_TRUE;

1233

} else {

1234

ca = c->a;

1235

ierr = PetscArrayzero(ca,ci[cm]+1)PetscMemzero(ca,(ci[cm]+1)*sizeof(*(ca)));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1235,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1236

}

1237

1238

if (abt->usecoloring) {

1239

MatTransposeColoring matcoloring = abt->matcoloring;

1240

Mat Bt_dense,C_dense = abt->ABt_den;

1241

1242

/* Get Bt_dense by Apply MatTransposeColoring to B */

1243

Bt_dense = abt->Bt_den;

1244

ierr = MatTransColoringApplySpToDen(matcoloring,B,Bt_dense);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1244,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1245

1246

/* C_dense = A*Bt_dense */

1247

ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,Bt_dense,C_dense);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1247,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1248

1249

/* Recover C from C_dense */

1250

ierr = MatTransColoringApplyDenToSp(matcoloring,C_dense,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1250,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1251

1252

}

1253

1254

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

1255

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

1256

acol = aj + ai[i];

1257

aval = aa + ai[i];

1258

cnzi = ci[i+1] - ci[i];

1259

ccol = cj + ci[i];

1260

cval = ca + ci[i];

1261

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

1262

brow = ccol[j];

1263

bnzj = bi[brow+1] - bi[brow];

1264

bcol = bj + bi[brow];

1265

bval = ba + bi[brow];

1266

1267

/* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */

1268

nexta = 0; nextb = 0;

1269

while (nexta<anzi && nextb<bnzj) {

1270

while (nexta < anzi && acol[nexta] < bcol[nextb]) nexta++;

1271

if (nexta == anzi) break;

1272

while (nextb < bnzj && acol[nexta] > bcol[nextb]) nextb++;

1273

if (nextb == bnzj) break;

1274

if (acol[nexta] == bcol[nextb]) {

1275

cval[j] += aval[nexta]*bval[nextb];

1276

nexta++; nextb++;

1277

flops += 2;

1278

}

1279

}

1280

}

1281

}

1282

ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1282,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1283

ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1283,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1284

ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1284,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1285

1286

}

1287

1288

PetscErrorCode MatDestroy_SeqAIJ_MatTransMatMult(Mat A)

1289

{

1290

PetscErrorCode ierr;

1291

Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;

1292

Mat_MatTransMatMult *atb = a->atb;

1293

1294

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1294; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1295

if (atb) {

1296

ierr = MatDestroy(&atb->At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1296,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1297

ierr = (*atb->destroy)(A);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1297,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1298

}

1299

ierr = PetscFree(atb)((*PetscTrFree)((void*)(atb),1299,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((atb) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1299,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1300

1301

}

1302

1303

PetscErrorCode MatTransposeMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)

1304

{

1305

PetscErrorCode ierr;

1306

const char *algTypes[2] = {"matmatmult","outerproduct"};

1307

PetscInt alg=0; /* set default algorithm */

1308

Mat At;

1309

Mat_MatTransMatMult *atb;

1310

Mat_SeqAIJ *c;

1311

1312

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1312; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1313

if (scall == MAT_INITIAL_MATRIX) {

1314

ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatTransposeMatMult","Mat")0; do { PetscOptionItems PetscOptionsObjectBase; PetscOptionItems
*PetscOptionsObject = &PetscOptionsObjectBase; PetscMemzero
(PetscOptionsObject,sizeof(PetscOptionItems)); for (PetscOptionsObject
->count=(PetscOptionsPublish?-1:1); PetscOptionsObject->
count<2; PetscOptionsObject->count++) { PetscErrorCode _5_ierr
= PetscOptionsBegin_Private(PetscOptionsObject,PetscObjectComm
((PetscObject)A),((PetscObject)A)->prefix,"MatTransposeMatMult"
,"Mat");do {if (__builtin_expect(!!(_5_ierr),0)) return PetscError
(((MPI_Comm)0x44000001),1314,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1314,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1315

ierr = PetscOptionsEList("-mattransposematmult_via","Algorithmic approach","MatTransposeMatMult",algTypes,2,algTypes[0],&alg,NULL)PetscOptionsEList_Private(PetscOptionsObject,"-mattransposematmult_via"
,"Algorithmic approach","MatTransposeMatMult",algTypes,2,algTypes
[0],&alg,((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1315,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1316

ierr = PetscOptionsEnd()_5_ierr = PetscOptionsEnd_Private(PetscOptionsObject);do {if (
__builtin_expect(!!(_5_ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1316,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);}} while (0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1316,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1317

1318

switch (alg) {

1319

case 1:

1320

ierr = MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1320,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1321

break;

1322

default:

1323

ierr = PetscNew(&atb)PetscMallocA(1,PETSC_TRUE,1323,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(1)*sizeof(**((&atb))),((&atb)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1323,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1324

ierr = MatTranspose_SeqAIJ(A,MAT_INITIAL_MATRIX,&At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1324,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1325

ierr = MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_INITIAL_MATRIX,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1325,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1326

1327

c = (Mat_SeqAIJ*)(*C)->data;

1328

c->atb = atb;

1329

atb->At = At;

1330

atb->destroy = (*C)->ops->destroy;

1331

(*C)->ops->destroy = MatDestroy_SeqAIJ_MatTransMatMult;

1332

1333

break;

1334

}

1335

}

1336

if (alg) {

1337

ierr = (*(*C)->ops->mattransposemultnumeric)(A,B,*C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1337,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1338

} else if (!alg && scall == MAT_REUSE_MATRIX) {

1339

c = (Mat_SeqAIJ*)(*C)->data;

1340

atb = c->atb;

1341

At = atb->At;

1342

ierr = MatTranspose_SeqAIJ(A,MAT_REUSE_MATRIX,&At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1342,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1343

ierr = MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_REUSE_MATRIX,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1343,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1344

}

1345

1346

}

1347

1348

PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)

1349

{

1350

PetscErrorCode ierr;

1351

Mat At;

1352

PetscInt *ati,*atj;

1353

1354

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1354; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1355

/* create symbolic At */

1356

ierr = MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1356,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1357

ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF((MPI_Comm)0x44000001),A->cmap->n,A->rmap->n,ati,atj,NULL((void*)0),&At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1357,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1358

ierr = MatSetBlockSizes(At,PetscAbs(A->cmap->bs)(((A->cmap->bs) >= 0) ? (A->cmap->bs) : (-(A->
cmap->bs))),PetscAbs(B->cmap->bs)(((B->cmap->bs) >= 0) ? (B->cmap->bs) : (-(B->
cmap->bs))));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1358,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1359

ierr = MatSetType(At,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1359,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1360

1361

/* get symbolic C=At*B */

1362

ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1362,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1363

1364

/* clean up */

1365

ierr = MatDestroy(&At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1365,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1366

ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1366,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1367

1368

(*C)->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ;

1369

1370

}

1371

1372

PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)

1373

{

1374

PetscErrorCode ierr;

1375

Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;

1376

PetscInt am =A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb;

1377

PetscInt cm =C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k;

1378

PetscLogDouble flops=0.0;

1379

MatScalar *aa =a->a,*ba,*ca,*caj;

1380

1381

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1381; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1382

if (!c->a) {

1383

ierr = PetscCalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_TRUE,1383,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[cm]+1)*sizeof(**(&ca)),(&ca));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1383,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1384

1385

c->a = ca;

1386

c->free_a = PETSC_TRUE;

1387

} else {

1388

ca = c->a;

1389

ierr = PetscArrayzero(ca,ci[cm])PetscMemzero(ca,(ci[cm])*sizeof(*(ca)));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1389,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1390

}

1391

1392

/* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */

1393

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

1394

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

1395

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

1396

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

1397

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

1398

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

1399

nextb = 0;

1400

crow = *aj++;

1401

cjj = cj + ci[crow];

1402

caj = ca + ci[crow];

1403

/* perform sparse axpy operation. Note cjj includes bj. */

1404

for (k=0; nextb<bnzi; k++) {

1405

if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */

1406

caj[k] += (*aa)*(*(ba+nextb));

1407

nextb++;

1408

}

1409

}

1410

flops += 2*bnzi;

1411

aa++;

1412

}

1413

}

1414

1415

/* Assemble the final matrix and clean up */

1416

ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1416,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1417

ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1417,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1418

ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1418,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1419

1420

}

1421

1422

PETSC_INTERNextern __attribute__((visibility ("hidden"))) PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)

1423

{

1424

PetscErrorCode ierr;

1425

1426

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1426; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1427

if (scall == MAT_INITIAL_MATRIX) {

1428

ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultSymbolic].active) ? (*PetscLogPLB)((MAT_MatMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1428,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1429

ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1429,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1430

ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultSymbolic].active) ? (*PetscLogPLE)((MAT_MatMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1430,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1431

}

1432

ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultNumeric].active) ? (*PetscLogPLB)((MAT_MatMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1432,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1433

ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1433,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1434

ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultNumeric].active) ? (*PetscLogPLE)((MAT_MatMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1434,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1435

1436

}

1437

1438

PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)

1439

{

1440

PetscErrorCode ierr;

1441

1442

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1442; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1443

ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1443,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1444

1445

(*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;

1446

1447

}

1448

1449

PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)

1450

{

1451

Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;

1452

Mat_SeqDense *bd = (Mat_SeqDense*)B->data;

1453

PetscErrorCode ierr;

1454

PetscScalar *c,*b,r1,r2,r3,r4,*c1,*c2,*c3,*c4,aatmp;

1455

const PetscScalar *aa,*b1,*b2,*b3,*b4;

1456

const PetscInt *aj;

1457

PetscInt cm=C->rmap->n,cn=B->cmap->n,bm=bd->lda,am=A->rmap->n;

1458

PetscInt am4=4*am,bm4=4*bm,col,i,j,n,ajtmp;

1459

1460

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1460; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1461

if (!cm || !cn) 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);

1462

if (B->rmap->n != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in A %D not equal rows in B %D\n",A->cmap->n,B->rmap->n)return PetscError(((MPI_Comm)0x44000001),1462,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,60,PETSC_ERROR_INITIAL,"Number columns in A %D not equal rows in B %D\n"
,A->cmap->n,B->rmap->n);

1463

if (A->rmap->n != C->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows in C %D not equal rows in A %D\n",C->rmap->n,A->rmap->n)return PetscError(((MPI_Comm)0x44000001),1463,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,60,PETSC_ERROR_INITIAL,"Number rows in C %D not equal rows in A %D\n"
,C->rmap->n,A->rmap->n);

1464

if (B->cmap->n != C->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in B %D not equal columns in C %D\n",B->cmap->n,C->cmap->n)return PetscError(((MPI_Comm)0x44000001),1464,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,60,PETSC_ERROR_INITIAL,"Number columns in B %D not equal columns in C %D\n"
,B->cmap->n,C->cmap->n);

1465

b = bd->v;

1466

ierr = MatDenseGetArray(C,&c);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1466,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1467

b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;

1468

c1 = c; c2 = c1 + am; c3 = c2 + am; c4 = c3 + am;

1469

for (col=0; col<cn-4; col += 4) { /* over columns of C */

1470

for (i=0; i<am; i++) { /* over rows of C in those columns */

1471

r1 = r2 = r3 = r4 = 0.0;

1472

n = a->i[i+1] - a->i[i];

1473

aj = a->j + a->i[i];

1474

aa = a->a + a->i[i];

1475

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

1476

aatmp = aa[j]; ajtmp = aj[j];

1477

r1 += aatmp*b1[ajtmp];

1478

r2 += aatmp*b2[ajtmp];

1479

r3 += aatmp*b3[ajtmp];

1480

r4 += aatmp*b4[ajtmp];

1481

}

1482

c1[i] = r1;

1483

c2[i] = r2;

1484

c3[i] = r3;

1485

c4[i] = r4;

1486

}

1487

b1 += bm4; b2 += bm4; b3 += bm4; b4 += bm4;

1488

c1 += am4; c2 += am4; c3 += am4; c4 += am4;

1489

}

1490

for (; col<cn; col++) { /* over extra columns of C */

1491

for (i=0; i<am; i++) { /* over rows of C in those columns */

1492

r1 = 0.0;

1493

n = a->i[i+1] - a->i[i];

1494

aj = a->j + a->i[i];

1495

aa = a->a + a->i[i];

1496

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

1497

r1 += aa[j]*b1[aj[j]];

1498

}

1499

c1[i] = r1;

1500

}

1501

b1 += bm;

1502

c1 += am;

1503

}

1504

ierr = PetscLogFlops(cn*(2.0*a->nz));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1504,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1505

ierr = MatDenseRestoreArray(C,&c);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1505,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1506

ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1506,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1507

ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1507,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1508

1509

}

1510

1511

1512

Note very similar to MatMult_SeqAIJ(), should generate both codes from same base

1513

1514

PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)

1515

{

1516

Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;

1517

Mat_SeqDense *bd = (Mat_SeqDense*)B->data;

1518

PetscErrorCode ierr;

1519

PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4;

1520

MatScalar *aa;

1521

PetscInt cm = C->rmap->n, cn=B->cmap->n, bm=bd->lda, col, i,j,n,*aj, am = A->rmap->n,*ii,arm;

1522

PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam,*ridx;

1523

1524

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1524; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1525

1526

b = bd->v;

1527

ierr = MatDenseGetArray(C,&c);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1527,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1528

b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;

1529

1530

if (a->compressedrow.use) { /* use compressed row format */

1531

for (col=0; col<cn-4; col += 4) { /* over columns of C */

1532

colam = col*am;

1533

arm = a->compressedrow.nrows;

1534

ii = a->compressedrow.i;

1535

ridx = a->compressedrow.rindex;

1536

for (i=0; i<arm; i++) { /* over rows of C in those columns */

1537

r1 = r2 = r3 = r4 = 0.0;

1538

n = ii[i+1] - ii[i];

1539

aj = a->j + ii[i];

1540

aa = a->a + ii[i];

1541

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

1542

r1 += (*aa)*b1[*aj];

1543

r2 += (*aa)*b2[*aj];

1544

r3 += (*aa)*b3[*aj];

1545

r4 += (*aa++)*b4[*aj++];

1546

}

1547

c[colam + ridx[i]] += r1;

1548

c[colam + am + ridx[i]] += r2;

1549

c[colam + am2 + ridx[i]] += r3;

1550

c[colam + am3 + ridx[i]] += r4;

1551

}

1552

b1 += bm4;

1553

b2 += bm4;

1554

b3 += bm4;

1555

b4 += bm4;

1556

}

1557

for (; col<cn; col++) { /* over extra columns of C */

1558

colam = col*am;

1559

arm = a->compressedrow.nrows;

1560

ii = a->compressedrow.i;

1561

ridx = a->compressedrow.rindex;

1562

for (i=0; i<arm; i++) { /* over rows of C in those columns */

1563

r1 = 0.0;

1564

n = ii[i+1] - ii[i];

1565

aj = a->j + ii[i];

1566

aa = a->a + ii[i];

1567

1568

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

1569

r1 += (*aa++)*b1[*aj++];

1570

}

1571

c[colam + ridx[i]] += r1;

1572

}

1573

b1 += bm;

1574

}

1575

} else {

1576

for (col=0; col<cn-4; col += 4) { /* over columns of C */

1577

colam = col*am;

1578

for (i=0; i<am; i++) { /* over rows of C in those columns */

1579

r1 = r2 = r3 = r4 = 0.0;

1580

n = a->i[i+1] - a->i[i];

1581

aj = a->j + a->i[i];

1582

aa = a->a + a->i[i];

1583

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

1584

r1 += (*aa)*b1[*aj];

1585

r2 += (*aa)*b2[*aj];

1586

r3 += (*aa)*b3[*aj];

1587

r4 += (*aa++)*b4[*aj++];

1588

}

1589

c[colam + i] += r1;

1590

c[colam + am + i] += r2;

1591

c[colam + am2 + i] += r3;

1592

c[colam + am3 + i] += r4;

1593

}

1594

b1 += bm4;

1595

b2 += bm4;

1596

b3 += bm4;

1597

b4 += bm4;

1598

}

1599

for (; col<cn; col++) { /* over extra columns of C */

1600

colam = col*am;

1601

for (i=0; i<am; i++) { /* over rows of C in those columns */

1602

r1 = 0.0;

1603

n = a->i[i+1] - a->i[i];

1604

aj = a->j + a->i[i];

1605

aa = a->a + a->i[i];

1606

1607

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

1608

r1 += (*aa++)*b1[*aj++];

1609

}

1610

c[colam + i] += r1;

1611

}

1612

b1 += bm;

1613

}

1614

}

1615

ierr = PetscLogFlops(cn*2.0*a->nz);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1615,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1616

ierr = MatDenseRestoreArray(C,&c);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1616,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1617

1618

}

1619

1620

PetscErrorCode MatTransColoringApplySpToDen_SeqAIJ(MatTransposeColoring coloring,Mat B,Mat Btdense)

1621

{

1622

PetscErrorCode ierr;

1623

Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;

1624

Mat_SeqDense *btdense = (Mat_SeqDense*)Btdense->data;

1625

PetscInt *bi = b->i,*bj=b->j;

1626

PetscInt m = Btdense->rmap->n,n=Btdense->cmap->n,j,k,l,col,anz,*btcol,brow,ncolumns;

1627

MatScalar *btval,*btval_den,*ba=b->a;

1628

PetscInt *columns=coloring->columns,*colorforcol=coloring->colorforcol,ncolors=coloring->ncolors;

1629

1630

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1630; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1631

btval_den=btdense->v;

1632

ierr = PetscArrayzero(btval_den,m*n)PetscMemzero(btval_den,(m*n)*sizeof(*(btval_den)));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1632,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1633

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

1634

ncolumns = coloring->ncolumns[k];

1635

for (l=0; l<ncolumns; l++) { /* insert a row of B to a column of Btdense */

1636

col = *(columns + colorforcol[k] + l);

1637

btcol = bj + bi[col];

1638

btval = ba + bi[col];

1639

anz = bi[col+1] - bi[col];

1640

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

1641

brow = btcol[j];

1642

btval_den[brow] = btval[j];

1643

}

1644

}

1645

btval_den += m;

1646

}

1647

1648

}

1649

1650

PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring,Mat Cden,Mat Csp)

1651

{

1652

PetscErrorCode ierr;

1653

Mat_SeqAIJ *csp = (Mat_SeqAIJ*)Csp->data;

1654

const PetscScalar *ca_den,*ca_den_ptr;

1655

PetscScalar *ca=csp->a;

1656

PetscInt k,l,m=Cden->rmap->n,ncolors=matcoloring->ncolors;

1657

PetscInt brows=matcoloring->brows,*den2sp=matcoloring->den2sp;

1658

PetscInt nrows,*row,*idx;

1659

PetscInt *rows=matcoloring->rows,*colorforrow=matcoloring->colorforrow;

1660

1661

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1661; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1662

ierr = MatDenseGetArrayRead(Cden,&ca_den);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1662,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1663

1664

if (brows > 0) {

1665

PetscInt *lstart,row_end,row_start;

1666

lstart = matcoloring->lstart;

1667

ierr = PetscArrayzero(lstart,ncolors)PetscMemzero(lstart,(ncolors)*sizeof(*(lstart)));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1667,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1668

1669

row_end = brows;

1670

if (row_end > m) row_end = m;

1671

for (row_start=0; row_start<m; row_start+=brows) { /* loop over row blocks of Csp */

1672

ca_den_ptr = ca_den;

1673

for (k=0; k<ncolors; k++) { /* loop over colors (columns of Cden) */

1674

nrows = matcoloring->nrows[k];

1675

row = rows + colorforrow[k];

1676

idx = den2sp + colorforrow[k];

1677

for (l=lstart[k]; l<nrows; l++) {

1678

if (row[l] >= row_end) {

1679

lstart[k] = l;

1680

break;

1681

} else {

1682

ca[idx[l]] = ca_den_ptr[row[l]];

1683

}

1684

}

1685

ca_den_ptr += m;

1686

}

1687

row_end += brows;

1688

if (row_end > m) row_end = m;

1689

}

1690

} else { /* non-blocked impl: loop over columns of Csp - slow if Csp is large */

1691

ca_den_ptr = ca_den;

1692

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

1693

nrows = matcoloring->nrows[k];

1694

row = rows + colorforrow[k];

1695

idx = den2sp + colorforrow[k];

1696

for (l=0; l<nrows; l++) {

1697

ca[idx[l]] = ca_den_ptr[row[l]];

1698

}

1699

ca_den_ptr += m;

1700

}

1701

}

1702

1703

ierr = MatDenseRestoreArrayRead(Cden,&ca_den);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1703,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1704

#if defined(PETSC_USE_INFO1)

1705

if (matcoloring->brows > 0) {

1706

ierr = PetscInfo1(Csp,"Loop over %D row blocks for den2sp\n",brows)PetscInfo_Private(__func__,Csp,"Loop over %D row blocks for den2sp\n"
,brows);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1706,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1707

} else {

1708

ierr = PetscInfo(Csp,"Loop over colors/columns of Cden, inefficient for large sparse matrix product \n")PetscInfo_Private(__func__,Csp,"Loop over colors/columns of Cden, inefficient for large sparse matrix product \n"
);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1708,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1709

}

1710

#endif

1711

1712

}

1713

1714

PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatTransposeColoring c)

1715

{

1716

PetscErrorCode ierr;

1717

PetscInt i,n,nrows,Nbs,j,k,m,ncols,col,cm;

1718

const PetscInt *is,*ci,*cj,*row_idx;

1719

PetscInt nis = iscoloring->n,*rowhit,bs = 1;

1720

IS *isa;

1721

Mat_SeqAIJ *csp = (Mat_SeqAIJ*)mat->data;

1722

PetscInt *colorforrow,*rows,*rows_i,*idxhit,*spidx,*den2sp,*den2sp_i;

1723

PetscInt *colorforcol,*columns,*columns_i,brows;

1724

PetscBool flg;

1725

1726

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1726; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1727

ierr = ISColoringGetIS(iscoloring,PETSC_USE_POINTER,PETSC_IGNORE((void*)0),&isa);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1727,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1728

1729

/* bs >1 is not being tested yet! */

1730

Nbs = mat->cmap->N/bs;

1731

c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */

1732

c->N = Nbs;

1733

c->m = c->M;

1734

c->rstart = 0;

1735

c->brows = 100;

1736

1737

c->ncolors = nis;

1738

ierr = PetscMalloc3(nis,&c->ncolumns,nis,&c->nrows,nis+1,&colorforrow)PetscMallocA(3,PETSC_FALSE,1738,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(nis)*sizeof(**(&c->ncolumns)),(&c->ncolumns
),(size_t)(nis)*sizeof(**(&c->nrows)),(&c->nrows
),(size_t)(nis+1)*sizeof(**(&colorforrow)),(&colorforrow
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1738,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1739

ierr = PetscMalloc1(csp->nz+1,&rows)PetscMallocA(1,PETSC_FALSE,1739,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(csp->nz+1)*sizeof(**(&rows)),(&rows));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1739,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1740

ierr = PetscMalloc1(csp->nz+1,&den2sp)PetscMallocA(1,PETSC_FALSE,1740,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(csp->nz+1)*sizeof(**(&den2sp)),(&den2sp)
);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1740,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1741

1742

brows = c->brows;

1743

ierr = PetscOptionsGetInt(NULL((void*)0),NULL((void*)0),"-matden2sp_brows",&brows,&flg);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1743,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1744

if (flg) c->brows = brows;

1745

if (brows > 0) {

1746

ierr = PetscMalloc1(nis+1,&c->lstart)PetscMallocA(1,PETSC_FALSE,1746,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(nis+1)*sizeof(**(&c->lstart)),(&c->lstart
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1746,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1747

}

1748

1749

colorforrow[0] = 0;

1750

rows_i = rows;

1751

den2sp_i = den2sp;

1752

1753

ierr = PetscMalloc1(nis+1,&colorforcol)PetscMallocA(1,PETSC_FALSE,1753,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(nis+1)*sizeof(**(&colorforcol)),(&colorforcol
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1753,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1754

ierr = PetscMalloc1(Nbs+1,&columns)PetscMallocA(1,PETSC_FALSE,1754,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(Nbs+1)*sizeof(**(&columns)),(&columns));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1754,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1755

1756

colorforcol[0] = 0;

1757

columns_i = columns;

1758

1759

/* get column-wise storage of mat */

1760

ierr = MatGetColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1760,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1761

1762

cm = c->m;

1763

ierr = PetscMalloc1(cm+1,&rowhit)PetscMallocA(1,PETSC_FALSE,1763,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(cm+1)*sizeof(**(&rowhit)),(&rowhit));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1763,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1764

ierr = PetscMalloc1(cm+1,&idxhit)PetscMallocA(1,PETSC_FALSE,1764,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(cm+1)*sizeof(**(&idxhit)),(&idxhit));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1764,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1765

for (i=0; i<nis; i++) { /* loop over color */

1766

ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1766,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1767

ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1767,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1768

1769

c->ncolumns[i] = n;

1770

if (n) {

1771

ierr = PetscArraycpy(columns_i,is,n)((sizeof(*(columns_i)) != sizeof(*(is))) || PetscMemcpy(columns_i
,is,(n)*sizeof(*(columns_i))));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1771,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1772

}

1773

colorforcol[i+1] = colorforcol[i] + n;

1774

columns_i += n;

1775

1776

/* fast, crude version requires O(N*N) work */

1777

ierr = PetscArrayzero(rowhit,cm)PetscMemzero(rowhit,(cm)*sizeof(*(rowhit)));;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1777,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1778

1779

for (j=0; j<n; j++) { /* loop over columns*/

1780

col = is[j];

1781

row_idx = cj + ci[col];

1782

m = ci[col+1] - ci[col];

1783

for (k=0; k<m; k++) { /* loop over columns marking them in rowhit */

1784

idxhit[*row_idx] = spidx[ci[col] + k];

1785

rowhit[*row_idx++] = col + 1;

1786

}

1787

}

1788

/* count the number of hits */

1789

nrows = 0;

1790

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

1791

if (rowhit[j]) nrows++;

1792

}

1793

c->nrows[i] = nrows;

1794

colorforrow[i+1] = colorforrow[i] + nrows;

1795

1796

nrows = 0;

1797

for (j=0; j<cm; j++) { /* loop over rows */

1798

if (rowhit[j]) {

1799

rows_i[nrows] = j;

1800

den2sp_i[nrows] = idxhit[j];

1801

nrows++;

1802

}

1803

}

1804

den2sp_i += nrows;

1805

1806

ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1806,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1807

rows_i += nrows;

1808

}

1809

ierr = MatRestoreColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1809,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1810

ierr = PetscFree(rowhit)((*PetscTrFree)((void*)(rowhit),1810,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((rowhit) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1810,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1811

ierr = ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&isa);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1811,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1812

if (csp->nz != colorforrow[nis]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"csp->nz %d != colorforrow[nis] %d",csp->nz,colorforrow[nis])return PetscError(((MPI_Comm)0x44000001),1812,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,77,PETSC_ERROR_INITIAL,"csp->nz %d != colorforrow[nis] %d"
,csp->nz,colorforrow[nis]);

1813

1814

c->colorforrow = colorforrow;

1815

c->rows = rows;

1816

c->den2sp = den2sp;

1817

c->colorforcol = colorforcol;

1818

c->columns = columns;

1819

1820

ierr = PetscFree(idxhit)((*PetscTrFree)((void*)(idxhit),1820,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((idxhit) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1820,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1821

1822

}

1823

1824

/* The combine method has done the symbolic and numeric in the first phase, and so we just return here */

1825

PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Combined(Mat A,Mat B,Mat C)

1826

{

1827

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1827; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1828

/* Empty function */

1829

1830

}

1831

1832

/* This algorithm combines the symbolic and numeric phase of matrix-matrix multiplication. */

1833

PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ_Combined(Mat A,Mat B,PetscReal fill,Mat *C)

1834

{

1835

PetscErrorCode ierr;

1836

PetscLogDouble flops=0.0;

1837

Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;

1838

const PetscInt *ai=a->i,*bi=b->i;

1839

PetscInt *ci,*cj,*cj_i;

1840

PetscScalar *ca,*ca_i;

1841

PetscInt b_maxmemrow,c_maxmem,a_col;

1842

PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;

1843

PetscInt i,k,ndouble=0;

1844

PetscReal afill;

1845

PetscScalar *c_row_val_dense;

1846

PetscBool *c_row_idx_flags;

1847

PetscInt *aj_i=a->j;

1848

PetscScalar *aa_i=a->a;

1849

1850

PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1850; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1851

1852

/* Step 1: Determine upper bounds on memory for C and allocate memory */

1853

/* This should be enough for almost all matrices. If still more memory is needed, it is reallocated later. */

1854

c_maxmem = 8*(ai[am]+bi[bm]);

1855

b_maxmemrow = PetscMin(bi[bm],bn)(((bi[bm])<(bn)) ? (bi[bm]) : (bn));

1856

ierr = PetscMalloc1(am+1,&ci)PetscMallocA(1,PETSC_FALSE,1856,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+1)*sizeof(**(&ci)),(&ci));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1856,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1857

ierr = PetscCalloc1(bn,&c_row_val_dense)PetscMallocA(1,PETSC_TRUE,1857,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(bn)*sizeof(**(&c_row_val_dense)),(&c_row_val_dense
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1857,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1858

ierr = PetscCalloc1(bn,&c_row_idx_flags)PetscMallocA(1,PETSC_TRUE,1858,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(bn)*sizeof(**(&c_row_idx_flags)),(&c_row_idx_flags
));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1858,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1859

ierr = PetscMalloc1(c_maxmem,&cj)PetscMallocA(1,PETSC_FALSE,1859,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(c_maxmem)*sizeof(**(&cj)),(&cj));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1859,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1860

ierr = PetscMalloc1(c_maxmem,&ca)PetscMallocA(1,PETSC_FALSE,1860,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(c_maxmem)*sizeof(**(&ca)),(&ca));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1860,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1861

ca_i = ca;

1862

cj_i = cj;

1863

ci[0] = 0;

1864

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

1865

/* Step 2: Initialize the dense row vector for C */

1866

const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */

1867

PetscInt cnzi = 0;

1868

PetscInt *bj_i;

1869

PetscScalar *ba_i;

1870

/* If the number of indices in C so far + the max number of columns in the next row > c_maxmem -> allocate more memory

1871

Usually, there is enough memory in the first place, so this is not executed. */

1872

while (ci[i] + b_maxmemrow > c_maxmem) {

1873

c_maxmem *= 2;

1874

ndouble++;

1875

ierr = PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj)((*PetscTrRealloc)((sizeof(PetscInt)*c_maxmem),1875,__func__,
"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(void**)(&cj)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1875,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1876

ierr = PetscRealloc(sizeof(PetscScalar)*c_maxmem,&ca)((*PetscTrRealloc)((sizeof(PetscScalar)*c_maxmem),1876,__func__
,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,(void**)(&ca)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1876,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1877

}

1878

1879

/* Step 3: Do the numerical calculations */

1880

for (a_col=0; a_col<anzi; a_col++) { /* iterate over all non zero values in a row of A */

1881

PetscInt a_col_index = aj_i[a_col];

1882

const PetscInt bnzi = bi[a_col_index+1] - bi[a_col_index];

1883

flops += 2*bnzi;

1884

bj_i = b->j + bi[a_col_index]; /* points to the current row in bj */

1885

ba_i = b->a + bi[a_col_index]; /* points to the current row in ba */

1886

for (k=0; k<bnzi; ++k) { /* iterate over all non zeros of this row in B */

1887

if (c_row_idx_flags[bj_i[k]] == PETSC_FALSE) {

1888

cj_i[cnzi++] = bj_i[k];

1889

c_row_idx_flags[bj_i[k]] = PETSC_TRUE;

1890

}

1891

c_row_val_dense[bj_i[k]] += aa_i[a_col] * ba_i[k];

1892

}

1893

}

1894

1895

/* Sort array */

1896

ierr = PetscSortInt(cnzi,cj_i);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1896,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1897

/* Step 4 */

1898

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

1899

ca_i[k] = c_row_val_dense[cj_i[k]];

1900

c_row_val_dense[cj_i[k]] = 0.;

1901

c_row_idx_flags[cj_i[k]] = PETSC_FALSE;

1902

}

1903

/* terminate current row */

1904

aa_i += anzi;

1905

aj_i += anzi;

1906

ca_i += cnzi;

1907

cj_i += cnzi;

1908

ci[i+1] = ci[i] + cnzi;

1909

flops += cnzi;

1910

}

1911

1912

/* Step 5 */

1913

/* Create the new matrix */

1914

ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1914,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1915

ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1915,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1916

ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1916,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1917

1918

/* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */

1919

/* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */

1920

c = (Mat_SeqAIJ*)((*C)->data);

1921

c->a = ca;

1922

c->free_a = PETSC_TRUE;

1923

c->free_ij = PETSC_TRUE;

1924

c->nonew = 0;

1925

1926

/* set MatInfo */

1927

afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;

1928

if (afill < 1.0) afill = 1.0;

1929

c->maxnz = ci[am];

1930

c->nz = ci[am];

1931

(*C)->info.mallocs = ndouble;

1932

(*C)->info.fill_ratio_given = fill;

1933

(*C)->info.fill_ratio_needed = afill;

1934

(*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Combined;

1935

1936

ierr = PetscFree(c_row_val_dense)((*PetscTrFree)((void*)(c_row_val_dense),1936,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((c_row_val_dense) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1936,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1937

ierr = PetscFree(c_row_idx_flags)((*PetscTrFree)((void*)(c_row_idx_flags),1937,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
) || ((c_row_idx_flags) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1937,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1938

1939

ierr = MatAssemblyBegin(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1939,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1940

ierr = MatAssemblyEnd(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1940,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1941

ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1941,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1942

1943

}