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

Bug Summary

File:	mat/impls/aij/seq/matmatmult.c
Warning:	line 685, 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),24,__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),25,__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),26,__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] = 34; 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),36,__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),38,__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] = 55; 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),56,__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),56,__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),57,__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),58,__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),58,__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),61,__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),64,__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),67,__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),70,__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),73,__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),76,__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),79,__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),83,__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),87,__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;

PetscReal afill;

100

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

101

PetscTable ta;

102

PetscBT lnkbt;

103

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

104

105

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] = 105; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

106

/* Get ci and cj */

107

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

108

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

109

/* free space for accumulating nonzero column info */

110

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,110,__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),110,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

111

ci[0] = 0;

112

113

/* create and initialize a linked list */

114

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

115

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

116

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

117

ierr = PetscTableDestroy(&ta);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

119

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

120

121

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

122

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

123

124

current_space = free_space;

125

126

/* Determine ci and cj */

127

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

128

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

129

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

130

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

131

brow = aj[j];

132

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

133

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

134

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

135

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

136

}

137

cnzi = lnk[0];

138

139

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

140

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

141

if (current_space->local_remaining<cnzi) {

142

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

143

ndouble++;

144

}

145

146

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

147

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

148

149

current_space->array += cnzi;

150

current_space->local_used += cnzi;

151

current_space->local_remaining -= cnzi;

152

153

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

154

}

155

156

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

157

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

158

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

159

ierr = PetscMalloc1(ci[am]+1,&cj)PetscMallocA(1,PETSC_FALSE,159,__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),159,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

160

ierr = PetscFreeSpaceContiguous(&free_space,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 = PetscLLCondensedDestroy(lnk,lnkbt);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

163

/* put together the new symbolic matrix */

164

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),164,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

165

ierr = MatSetBlockSizesFromMats(*C,A,B);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 = MatSetType(*C,((PetscObject)A)->type_name);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

168

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

169

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

170

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

171

c->free_a = PETSC_FALSE;

172

c->free_ij = PETSC_TRUE;

173

c->nonew = 0;

174

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

175

176

/* set MatInfo */

177

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

178

if (afill < 1.0) afill = 1.0;

179

c->maxnz = ci[am];

180

c->nz = ci[am];

181

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

182

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

183

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

184

185

#if defined(PETSC_USE_INFO1)

186

if (ci[am]) {

187

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),187,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

188

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),188,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

189

} else {

190

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),190,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

191

}

192

#endif

193

194

}

195

196

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

197

{

198

PetscErrorCode ierr;

199

PetscLogDouble flops=0.0;

200

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

201

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

202

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

203

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

204

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

205

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

206

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

207

PetscScalar *ab_dense;

208

209

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] = 209; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

210

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

211

ierr = PetscMalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_FALSE,211,__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),211,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

212

c->a = ca;

213

c->free_a = PETSC_TRUE;

214

} else {

215

ca = c->a;

216

}

217

if (!c->matmult_abdense) {

218

ierr = PetscCalloc1(B->cmap->N,&ab_dense)PetscMallocA(1,PETSC_TRUE,218,__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),218,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

219

c->matmult_abdense = ab_dense;

220

} else {

221

ab_dense = c->matmult_abdense;

222

}

223

224

/* clean old values in C */

225

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

226

/* Traverse A row-wise. */

227

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

228

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

229

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

230

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

231

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

232

brow = aj[j];

233

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

234

bjj = bj + bi[brow];

235

baj = ba + bi[brow];

236

/* perform dense axpy */

237

valtmp = aa[j];

238

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

239

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

240

}

241

flops += 2*bnzi;

242

}

243

aj += anzi; aa += anzi;

244

245

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

246

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

247

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

248

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

249

}

250

flops += cnzi;

251

cj += cnzi; ca += cnzi;

252

}

253

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

254

ierr = MatAssemblyEnd(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 = PetscLogFlops(flops);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

257

}

258

259

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

260

{

261

PetscErrorCode ierr;

262

PetscLogDouble flops=0.0;

263

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

264

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

265

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

266

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

267

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

268

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

269

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

270

PetscInt nextb;

271

272

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] = 272; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

273

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

274

ierr = PetscMalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_FALSE,274,__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),274,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

275

c->a = ca;

276

c->free_a = PETSC_TRUE;

277

}

278

279

/* clean old values in C */

280

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

281

/* Traverse A row-wise. */

282

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

283

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

284

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

285

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

286

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

287

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

288

brow = aj[j];

289

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

290

bjj = bj + bi[brow];

291

baj = ba + bi[brow];

292

/* perform sparse axpy */

293

valtmp = aa[j];

294

nextb = 0;

295

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

296

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

297

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

298

}

299

}

300

flops += 2*bnzi;

301

}

302

aj += anzi; aa += anzi;

303

cj += cnzi; ca += cnzi;

304

}

305

306

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

307

ierr = MatAssemblyEnd(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 = PetscLogFlops(flops);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

310

}

311

312

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

313

{

314

PetscErrorCode ierr;

315

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

316

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

317

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

318

MatScalar *ca;

319

PetscReal afill;

320

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

321

PetscTable ta;

322

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

323

324

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] = 324; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

325

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

326

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

327

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

328

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,328,__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),328,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

329

ci[0] = 0;

330

331

/* create and initialize a linked list */

332

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

333

334

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

335

ierr = PetscTableDestroy(&ta);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

337

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

338

339

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

340

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

341

current_space = free_space;

342

343

/* Determine ci and cj */

344

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

345

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

346

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

347

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

348

brow = aj[j];

349

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

350

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

351

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

352

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

353

}

354

cnzi = lnk[1];

355

356

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

357

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

358

if (current_space->local_remaining<cnzi) {

359

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

360

ndouble++;

361

}

362

363

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

364

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

365

366

current_space->array += cnzi;

367

current_space->local_used += cnzi;

368

current_space->local_remaining -= cnzi;

369

370

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

371

}

372

373

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

374

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

375

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

376

ierr = PetscMalloc1(ci[am]+1,&cj)PetscMallocA(1,PETSC_FALSE,376,__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),376,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

377

ierr = PetscFreeSpaceContiguous(&free_space,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 = PetscLLCondensedDestroy_fast(lnk);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

380

/* Allocate space for ca */

381

ierr = PetscCalloc1(ci[am]+1,&ca)PetscMallocA(1,PETSC_TRUE,381,__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),381,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

382

383

/* put together the new symbolic matrix */

384

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

385

ierr = MatSetBlockSizesFromMats(*C,A,B);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 = MatSetType(*C,((PetscObject)A)->type_name);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

388

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

389

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

390

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

391

c->free_a = PETSC_TRUE;

392

c->free_ij = PETSC_TRUE;

393

c->nonew = 0;

394

395

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

396

397

/* set MatInfo */

398

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

399

if (afill < 1.0) afill = 1.0;

400

c->maxnz = ci[am];

401

c->nz = ci[am];

402

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

403

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

404

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

405

406

#if defined(PETSC_USE_INFO1)

407

if (ci[am]) {

408

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),408,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

409

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),409,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

410

} else {

411

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),411,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

412

}

413

#endif

414

415

}

416

417

418

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

419

{

420

PetscErrorCode ierr;

421

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

422

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

423

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

424

MatScalar *ca;

425

PetscReal afill;

426

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

427

PetscTable ta;

428

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

429

430

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] = 430; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

431

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

432

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

433

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

434

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,434,__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),434,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

435

ci[0] = 0;

436

437

/* create and initialize a linked list */

438

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

439

440

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

441

ierr = PetscTableDestroy(&ta);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 = PetscLLCondensedCreate_Scalable(Crmax,&lnk);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

444

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

445

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

446

current_space = free_space;

447

448

/* Determine ci and cj */

449

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

450

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

451

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

452

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

453

brow = aj[j];

454

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

455

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

456

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

457

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

458

}

459

cnzi = lnk[0];

460

461

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

462

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

463

if (current_space->local_remaining<cnzi) {

464

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

465

ndouble++;

466

}

467

468

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

469

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

470

471

current_space->array += cnzi;

472

current_space->local_used += cnzi;

473

current_space->local_remaining -= cnzi;

474

475

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

476

}

477

478

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

479

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

480

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

481

ierr = PetscMalloc1(ci[am]+1,&cj)PetscMallocA(1,PETSC_FALSE,481,__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),481,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

482

ierr = PetscFreeSpaceContiguous(&free_space,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 = PetscLLCondensedDestroy_Scalable(lnk);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

485

/* Allocate space for ca */

486

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

487

ierr = PetscCalloc1(ci[am]+1,&ca)PetscMallocA(1,PETSC_TRUE,487,__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),487,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

488

489

/* put together the new symbolic matrix */

490

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

491

ierr = MatSetBlockSizesFromMats(*C,A,B);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 = MatSetType(*C,((PetscObject)A)->type_name);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

494

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

495

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

496

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

497

c->free_a = PETSC_TRUE;

498

c->free_ij = PETSC_TRUE;

499

c->nonew = 0;

500

501

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

502

503

/* set MatInfo */

504

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

505

if (afill < 1.0) afill = 1.0;

506

c->maxnz = ci[am];

507

c->nz = ci[am];

508

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

509

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

510

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

511

512

#if defined(PETSC_USE_INFO1)

513

if (ci[am]) {

514

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),514,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

515

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),515,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

516

} else {

517

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),517,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

518

}

519

#endif

520

521

}

522

523

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

524

{

525

PetscErrorCode ierr;

526

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

527

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

528

PetscInt *ci,*cj,*bb;

529

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

530

PetscReal afill;

531

PetscInt i,j,col,ndouble = 0;

532

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

533

PetscHeap h;

534

535

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] = 535; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

536

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

537

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

538

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

539

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,539,__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),539,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

540

ci[0] = 0;

541

542

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

543

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

544

current_space = free_space;

545

546

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

547

ierr = PetscMalloc1(a->rmax,&bb)PetscMallocA(1,PETSC_FALSE,547,__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),547,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

548

549

/* Determine ci and cj */

550

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

551

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

552

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

553

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

554

/* Populate the min heap */

555

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

556

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

557

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

558

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

559

}

560

}

561

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

562

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

563

while (j >= 0) {

564

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

565

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),565,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

566

ndouble++;

567

}

568

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

569

current_space->local_used++;

570

current_space->local_remaining--;

571

ci[i+1]++;

572

573

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

574

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),574,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);}

575

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

576

PetscInt j2,col2;

577

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

578

if (col2 != col) break;

579

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

580

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),580,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);}

581

}

582

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

583

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

584

ierr = PetscHeapPop(h,&j,&col);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

}

586

}

587

ierr = PetscFree(bb)((*PetscTrFree)((void*)(bb),587,__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),587,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

588

ierr = PetscHeapDestroy(&h);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

590

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

591

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

592

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

593

ierr = PetscMalloc1(ci[am],&cj)PetscMallocA(1,PETSC_FALSE,593,__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),593,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

594

ierr = PetscFreeSpaceContiguous(&free_space,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

596

/* put together the new symbolic matrix */

597

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),597,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

598

ierr = MatSetBlockSizesFromMats(*C,A,B);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 = MatSetType(*C,((PetscObject)A)->type_name);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

601

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

602

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

603

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

604

c->free_a = PETSC_TRUE;

605

c->free_ij = PETSC_TRUE;

606

c->nonew = 0;

607

608

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

609

610

/* set MatInfo */

611

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

612

if (afill < 1.0) afill = 1.0;

613

c->maxnz = ci[am];

614

c->nz = ci[am];

615

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

616

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

617

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

618

619

#if defined(PETSC_USE_INFO1)

620

if (ci[am]) {

621

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),621,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

622

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),622,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

623

} else {

624

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),624,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

625

}

626

#endif

627

628

}

629

630

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

631

{

632

PetscErrorCode ierr;

633

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

634

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

635

PetscInt *ci,*cj,*bb;

636

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

637

PetscReal afill;

638

PetscInt i,j,col,ndouble = 0;

639

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

640

PetscHeap h;

641

PetscBT bt;

642

643

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] = 643; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

644

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

645

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

646

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

647

ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,647,__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),647,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

648

ci[0] = 0;

649

650

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

651

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

652

653

current_space = free_space;

Value assigned to 'current_space'

→

654

655

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

656

ierr = PetscMalloc1(a->rmax,&bb)PetscMallocA(1,PETSC_FALSE,656,__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),656,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

657

ierr = PetscBTCreate(bn,&bt);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

659

/* Determine ci and cj */

660

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

→

661

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

662

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

663

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

←

Assuming pointer value is null

→

664

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

665

/* Populate the min heap */

666

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

←

Assuming 'j' is >= 'anzi'

→

←

Loop condition is false. Execution continues on line 678

→

←

Assuming 'j' is >= 'anzi'

→

←

Loop condition is false. Execution continues on line 678

→

667

PetscInt brow = acol[j];

668

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

669

PetscInt bcol = bj[bb[j]];

670

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

671

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

672

bb[j]++;

673

break;

674

}

675

}

676

}

677

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

678

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

679

while (j >= 0) {

←

Assuming 'j' is < 0

→

←

Loop condition is false. Execution continues on line 701

Assuming 'j' is >= 0

Loop condition is true. Entering loop body

→

680

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

←

Assuming the condition is false

→

←

Taking false branch

→

681

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

682

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),682,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

683

ndouble++;

684

}

685

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

←

Dereference of null pointer

686

current_space->local_used++;

687

current_space->local_remaining--;

688

ci[i+1]++;

689

690

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

691

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

692

PetscInt bcol = bj[bb[j]];

693

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

694

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

695

bb[j]++;

696

break;

697

}

698

}

699

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

700

}

701

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

←

Assuming 'fptr' is null

→

←

Taking false branch

→

702

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

703

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

704

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

705

}

706

}

707

ierr = PetscFree(bb)((*PetscTrFree)((void*)(bb),707,__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),707,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

708

ierr = PetscHeapDestroy(&h);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 = PetscBTDestroy(&bt);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

711

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

712

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

713

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

714

ierr = PetscMalloc1(ci[am],&cj)PetscMallocA(1,PETSC_FALSE,714,__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),714,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

715

ierr = PetscFreeSpaceContiguous(&free_space,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

717

/* put together the new symbolic matrix */

718

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),718,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

719

ierr = MatSetBlockSizesFromMats(*C,A,B);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 = MatSetType(*C,((PetscObject)A)->type_name);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

722

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

723

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

724

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

725

c->free_a = PETSC_TRUE;

726

c->free_ij = PETSC_TRUE;

727

c->nonew = 0;

728

729

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

730

731

/* set MatInfo */

732

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

733

if (afill < 1.0) afill = 1.0;

734

c->maxnz = ci[am];

735

c->nz = ci[am];

736

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

737

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

738

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

739

740

#if defined(PETSC_USE_INFO1)

741

if (ci[am]) {

742

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),742,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

743

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),743,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

744

} else {

745

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),745,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

746

}

747

#endif

748

749

}

750

751

752

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

753

{

754

PetscErrorCode ierr;

755

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

756

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

757

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

758

PetscInt c_maxmem,a_maxrownnz=0,a_rownnz;

759

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

760

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

761

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

762

PetscInt window[8];

763

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

764

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

765

PetscReal afill;

766

PetscInt *workj_L1,*workj_L2,*workj_L3;

767

PetscInt L1_nnz,L2_nnz;

768

769

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

770

Because of the way virtual memory works,

771

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

772

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] = 772; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

773

ierr = PetscMalloc1(am+1,&ci)PetscMallocA(1,PETSC_FALSE,773,__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),773,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

774

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

775

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

776

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

777

a_rownnz = 0;

778

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

779

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

780

if (a_rownnz > bn) {

781

a_rownnz = bn;

782

break;

783

}

784

}

785

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

786

}

787

/* temporary work areas for merging rows */

788

ierr = PetscMalloc1(a_maxrownnz*8,&workj_L1)PetscMallocA(1,PETSC_FALSE,788,__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),788,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

789

ierr = PetscMalloc1(a_maxrownnz*8,&workj_L2)PetscMallocA(1,PETSC_FALSE,789,__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),789,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

790

ierr = PetscMalloc1(a_maxrownnz,&workj_L3)PetscMallocA(1,PETSC_FALSE,790,__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),790,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

791

792

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

793

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

794

/* Step 2: Populate pattern for C */

795

ierr = PetscMalloc1(c_maxmem,&cj)PetscMallocA(1,PETSC_FALSE,795,__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),795,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

796

797

ci_nnz = 0;

798

ci[0] = 0;

799

worki_L1[0] = 0;

800

worki_L2[0] = 0;

801

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

802

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

803

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

804

rowsleft = anzi;

805

inputcol_L1 = acol;

806

L2_nnz = 0;

807

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

808

worki_L2[1] = 0;

809

outputi_nnz = 0;

810

811

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

812

while (ci_nnz+a_maxrownnz > c_maxmem) {

813

c_maxmem *= 2;

814

ndouble++;

815

ierr = PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj)((*PetscTrRealloc)((sizeof(PetscInt)*c_maxmem),815,__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),815,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

816

}

817

818

while (rowsleft) {

819

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

820

L1_nrows = 0;

821

L1_nnz = 0;

822

inputcol = inputcol_L1;

823

inputi = bi;

824

inputj = bj;

825

826

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

827

This helps with compiler unrolling, improving performance substantially.

828

Input: inputj inputi inputcol bn

829

Output: outputj outputi_nnz */

830

#define MatMatMultSymbolic_RowMergeMacro(ANNZ) \

831

window_min = bn; \

832

outputi_nnz = 0; \

833

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

834

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

835

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

836

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

837

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

838

} \

839

while (window_min < bn) { \

840

outputj[outputi_nnz++] = window_min; \

841

/* advance front and compute new minimum */ \

842

old_window_min = window_min; \

843

window_min = bn; \

844

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

845

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

846

brow_ptr[k]++; \

847

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

848

} \

849

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

850

} \

851

}

852

853

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

854

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

855

while (L1_rowsleft) {

856

outputi_nnz = 0;

857

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

858

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

859

860

switch (L1_rowsleft) {

861

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

862

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

863

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

864

inputcol += L1_rowsleft;

865

rowsleft -= L1_rowsleft;

866

L1_rowsleft = 0;

867

break;

868

case 2: MatMatMultSymbolic_RowMergeMacro(2);

869

inputcol += L1_rowsleft;

870

rowsleft -= L1_rowsleft;

871

L1_rowsleft = 0;

872

break;

873

case 3: MatMatMultSymbolic_RowMergeMacro(3);

874

inputcol += L1_rowsleft;

875

rowsleft -= L1_rowsleft;

876

L1_rowsleft = 0;

877

break;

878

case 4: MatMatMultSymbolic_RowMergeMacro(4);

879

inputcol += L1_rowsleft;

880

rowsleft -= L1_rowsleft;

881

L1_rowsleft = 0;

882

break;

883

case 5: MatMatMultSymbolic_RowMergeMacro(5);

884

inputcol += L1_rowsleft;

885

rowsleft -= L1_rowsleft;

886

L1_rowsleft = 0;

887

break;

888

case 6: MatMatMultSymbolic_RowMergeMacro(6);

889

inputcol += L1_rowsleft;

890

rowsleft -= L1_rowsleft;

891

L1_rowsleft = 0;

892

break;

893

case 7: MatMatMultSymbolic_RowMergeMacro(7);

894

inputcol += L1_rowsleft;

895

rowsleft -= L1_rowsleft;

896

L1_rowsleft = 0;

897

break;

898

default: MatMatMultSymbolic_RowMergeMacro(8);

899

inputcol += 8;

900

rowsleft -= 8;

901

L1_rowsleft -= 8;

902

break;

903

}

904

inputcol_L1 = inputcol;

905

L1_nnz += outputi_nnz;

906

worki_L1[++L1_nrows] = L1_nnz;

907

}

908

909

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

910

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

911

if (anzi > 8) {

912

inputi = worki_L1;

913

inputj = workj_L1;

914

inputcol = workcol;

915

outputi_nnz = 0;

916

917

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

918

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

919

920

switch (L1_nrows) {

921

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

922

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

923

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

924

break;

925

case 2: MatMatMultSymbolic_RowMergeMacro(2); break;

926

case 3: MatMatMultSymbolic_RowMergeMacro(3); break;

927

case 4: MatMatMultSymbolic_RowMergeMacro(4); break;

928

case 5: MatMatMultSymbolic_RowMergeMacro(5); break;

929

case 6: MatMatMultSymbolic_RowMergeMacro(6); break;

930

case 7: MatMatMultSymbolic_RowMergeMacro(7); break;

931

case 8: MatMatMultSymbolic_RowMergeMacro(8); break;

932

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),932,__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!"
);

933

}

934

L2_nnz += outputi_nnz;

935

worki_L2[++L2_nrows] = L2_nnz;

936

937

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

938

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

939

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

940

inputi = worki_L2;

941

inputj = workj_L2;

942

inputcol = workcol;

943

outputi_nnz = 0;

944

if (rowsleft) outputj = workj_L3;

945

else outputj = cj + ci_nnz;

946

switch (L2_nrows) {

947

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

948

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

949

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

950

break;

951

case 2: MatMatMultSymbolic_RowMergeMacro(2); break;

952

case 3: MatMatMultSymbolic_RowMergeMacro(3); break;

953

case 4: MatMatMultSymbolic_RowMergeMacro(4); break;

954

case 5: MatMatMultSymbolic_RowMergeMacro(5); break;

955

case 6: MatMatMultSymbolic_RowMergeMacro(6); break;

956

case 7: MatMatMultSymbolic_RowMergeMacro(7); break;

957

case 8: MatMatMultSymbolic_RowMergeMacro(8); break;

958

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),958,__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!"
);

959

}

960

L2_nrows = 1;

961

L2_nnz = outputi_nnz;

962

worki_L2[1] = outputi_nnz;

963

/* Copy to workj_L2 */

964

if (rowsleft) {

965

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

966

}

967

}

968

}

969

} /* while (rowsleft) */

970

#undef MatMatMultSymbolic_RowMergeMacro

971

972

/* terminate current row */

973

ci_nnz += outputi_nnz;

974

ci[i+1] = ci_nnz;

975

}

976

977

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

978

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),978,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

979

ierr = MatSetBlockSizesFromMats(*C,A,B);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 = MatSetType(*C,((PetscObject)A)->type_name);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

982

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

983

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

984

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

985

c->free_a = PETSC_TRUE;

986

c->free_ij = PETSC_TRUE;

987

c->nonew = 0;

988

989

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

990

991

/* set MatInfo */

992

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

993

if (afill < 1.0) afill = 1.0;

994

c->maxnz = ci[am];

995

c->nz = ci[am];

996

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

997

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

998

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

999

1000

#if defined(PETSC_USE_INFO1)

1001

if (ci[am]) {

1002

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),1002,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1003

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),1003,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1004

} else {

1005

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),1005,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1006

}

1007

#endif

1008

1009

/* Step 4: Free temporary work areas */

1010

ierr = PetscFree(workj_L1)((*PetscTrFree)((void*)(workj_L1),1010,__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),1010,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1011

ierr = PetscFree(workj_L2)((*PetscTrFree)((void*)(workj_L2),1011,__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),1011,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1012

ierr = PetscFree(workj_L3)((*PetscTrFree)((void*)(workj_L3),1012,__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),1012,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1013

1014

}

1015

1016

/* concatenate unique entries and then sort */

1017

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

1018

{

1019

PetscErrorCode ierr;

1020

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

1021

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

1022

PetscInt *ci,*cj;

1023

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

1024

PetscReal afill;

1025

PetscInt i,j,ndouble = 0;

1026

PetscSegBuffer seg,segrow;

1027

char *seen;

1028

1029

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] = 1029; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1030

ierr = PetscMalloc1(am+1,&ci)PetscMallocA(1,PETSC_FALSE,1030,__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),1030,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1031

ci[0] = 0;

1032

1033

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

1034

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

1035

ierr = PetscSegBufferCreate(sizeof(PetscInt),100,&segrow);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 = PetscCalloc1(bn,&seen)PetscMallocA(1,PETSC_TRUE,1036,__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),1036,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1037

1038

/* Determine ci and cj */

1039

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

1040

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

1041

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

1042

PetscInt packlen = 0,*PETSC_RESTRICT__restrict crow;

1043

/* Pack segrow */

1044

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

1045

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

1046

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

1047

PetscInt bcol = bj[k];

1048

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

1049

PetscInt *PETSC_RESTRICT__restrict slot;

1050

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

1051

*slot = bcol;

1052

seen[bcol] = 1;

1053

packlen++;

1054

}

1055

}

1056

}

1057

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

1058

ierr = PetscSegBufferExtractTo(segrow,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 = PetscSortInt(packlen,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

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

1061

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

1062

}

1063

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

1064

ierr = PetscFree(seen)((*PetscTrFree)((void*)(seen),1064,__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),1064,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1065

1066

/* Column indices are in the segmented buffer */

1067

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

1068

ierr = PetscSegBufferDestroy(&seg);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

1070

/* put together the new symbolic matrix */

1071

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),1071,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1072

ierr = MatSetBlockSizesFromMats(*C,A,B);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 = MatSetType(*C,((PetscObject)A)->type_name);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

1075

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

1076

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

1077

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

1078

c->free_a = PETSC_TRUE;

1079

c->free_ij = PETSC_TRUE;

1080

c->nonew = 0;

1081

1082

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

1083

1084

/* set MatInfo */

1085

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

1086

if (afill < 1.0) afill = 1.0;

1087

c->maxnz = ci[am];

1088

c->nz = ci[am];

1089

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

1090

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

1091

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

1092

1093

#if defined(PETSC_USE_INFO1)

1094

if (ci[am]) {

1095

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),1095,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1096

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),1096,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1097

} else {

1098

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),1098,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1099

}

1100

#endif

1101

1102

}

1103

1104

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

1105

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

1106

{

1107

PetscErrorCode ierr;

1108

1109

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] = 1109; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1110

if (scall == MAT_INITIAL_MATRIX) {

1111

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),1111,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1112

ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);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 = 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),1113,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1114

}

1115

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),1115,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1116

ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);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 = 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),1117,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1118

1119

}

1120

1121

PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(Mat A)

1122

{

1123

PetscErrorCode ierr;

1124

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

1125

Mat_MatMatTransMult *abt=a->abt;

1126

1127

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] = 1127; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1128

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

1129

ierr = MatTransposeColoringDestroy(&abt->matcoloring);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 = MatDestroy(&abt->Bt_den);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->ABt_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 = PetscFree(abt)((*PetscTrFree)((void*)(abt),1132,__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),1132,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1133

1134

}

1135

1136

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

1137

{

1138

PetscErrorCode ierr;

1139

Mat Bt;

1140

PetscInt *bti,*btj;

1141

Mat_MatMatTransMult *abt;

1142

Mat_SeqAIJ *c;

1143

1144

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] = 1144; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1145

/* create symbolic Bt */

1146

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

1147

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),1147,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1148

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),1148,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1149

ierr = MatSetType(Bt,((PetscObject)A)->type_name);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

1151

/* get symbolic C=A*Bt */

1152

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

1153

1154

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

1155

ierr = PetscNew(&abt)PetscMallocA(1,PETSC_TRUE,1155,__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),1155,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1156

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

1157

c->abt = abt;

1158

1159

abt->usecoloring = PETSC_FALSE;

1160

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

1161

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

1162

1163

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),1163,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1164

if (abt->usecoloring) {

1165

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

1166

MatTransposeColoring matcoloring;

1167

MatColoring coloring;

1168

ISColoring iscoloring;

1169

Mat Bt_dense,C_dense;

1170

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

1171

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

1172

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),1172,__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'"
);

1173

1174

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

1175

ierr = MatColoringSetDistance(coloring,2);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 = MatColoringSetType(coloring,MATCOLORINGSL"sl");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 = MatColoringSetFromOptions(coloring);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 = MatColoringApply(coloring,&iscoloring);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 = MatColoringDestroy(&coloring);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 = MatTransposeColoringCreate(*C,iscoloring,&matcoloring);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

1182

abt->matcoloring = matcoloring;

1183

1184

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

1185

1186

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

1187

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

1188

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),1188,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1189

ierr = MatSetType(Bt_dense,MATSEQDENSE"seqdense");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 = MatSeqDenseSetPreallocation(Bt_dense,NULL((void*)0));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

1192

Bt_dense->assembled = PETSC_TRUE;

1193

abt->Bt_den = Bt_dense;

1194

1195

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

1196

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),1196,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1197

ierr = MatSetType(C_dense,MATSEQDENSE"seqdense");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 = MatSeqDenseSetPreallocation(C_dense,NULL((void*)0));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

1200

Bt_dense->assembled = PETSC_TRUE;

1201

abt->ABt_den = C_dense;

1202

1203

#if defined(PETSC_USE_INFO1)

1204

{

1205

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

1206

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),1206,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1207

}

1208

#endif

1209

}

1210

/* clean up */

1211

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

1212

ierr = MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);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

1214

}

1215

1216

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

1217

{

1218

PetscErrorCode ierr;

1219

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

1220

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

1221

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

1222

PetscLogDouble flops=0.0;

1223

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

1224

Mat_MatMatTransMult *abt = c->abt;

1225

1226

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] = 1226; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1227

/* clear old values in C */

1228

if (!c->a) {

1229

ierr = PetscCalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_TRUE,1229,__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),1229,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1230

c->a = ca;

1231

c->free_a = PETSC_TRUE;

1232

} else {

1233

ca = c->a;

1234

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),1234,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1235

}

1236

1237

if (abt->usecoloring) {

1238

MatTransposeColoring matcoloring = abt->matcoloring;

1239

Mat Bt_dense,C_dense = abt->ABt_den;

1240

1241

/* Get Bt_dense by Apply MatTransposeColoring to B */

1242

Bt_dense = abt->Bt_den;

1243

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

1244

1245

/* C_dense = A*Bt_dense */

1246

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

1247

1248

/* Recover C from C_dense */

1249

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

1250

1251

}

1252

1253

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

1254

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

1255

acol = aj + ai[i];

1256

aval = aa + ai[i];

1257

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

1258

ccol = cj + ci[i];

1259

cval = ca + ci[i];

1260

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

1261

brow = ccol[j];

1262

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

1263

bcol = bj + bi[brow];

1264

bval = ba + bi[brow];

1265

1266

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

1267

nexta = 0; nextb = 0;

1268

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

1269

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

1270

if (nexta == anzi) break;

1271

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

1272

if (nextb == bnzj) break;

1273

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

1274

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

1275

nexta++; nextb++;

1276

flops += 2;

1277

}

1278

}

1279

}

1280

}

1281

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

1282

ierr = MatAssemblyEnd(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 = PetscLogFlops(flops);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

1285

}

1286

1287

PetscErrorCode MatDestroy_SeqAIJ_MatTransMatMult(Mat A)

1288

{

1289

PetscErrorCode ierr;

1290

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

1291

Mat_MatTransMatMult *atb = a->atb;

1292

1293

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] = 1293; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1294

if (atb) {

1295

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

1296

ierr = (*atb->destroy)(A);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

}

1298

ierr = PetscFree(atb)((*PetscTrFree)((void*)(atb),1298,__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),1298,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1299

1300

}

1301

1302

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

1303

{

1304

PetscErrorCode ierr;

1305

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

1306

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

1307

Mat At;

1308

Mat_MatTransMatMult *atb;

1309

Mat_SeqAIJ *c;

1310

1311

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] = 1311; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1312

if (scall == MAT_INITIAL_MATRIX) {

1313

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),1313,__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),1313,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1314

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),1314,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1315

ierr = PetscOptionsEnd()_5_ierr = PetscOptionsEnd_Private(PetscOptionsObject);do {if (
__builtin_expect(!!(_5_ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1315,__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),1315,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1316

1317

switch (alg) {

1318

case 1:

1319

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

1320

break;

1321

default:

1322

ierr = PetscNew(&atb)PetscMallocA(1,PETSC_TRUE,1322,__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),1322,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1323

ierr = MatTranspose_SeqAIJ(A,MAT_INITIAL_MATRIX,&At);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 = MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_INITIAL_MATRIX,fill,C);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

1326

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

1327

c->atb = atb;

1328

atb->At = At;

1329

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

1330

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

1331

1332

break;

1333

}

1334

}

1335

if (alg) {

1336

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

1337

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

1338

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

1339

atb = c->atb;

1340

At = atb->At;

1341

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

1342

ierr = MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_REUSE_MATRIX,fill,C);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

}

1344

1345

}

1346

1347

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

1348

{

1349

PetscErrorCode ierr;

1350

Mat At;

1351

PetscInt *ati,*atj;

1352

1353

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] = 1353; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1354

/* create symbolic At */

1355

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

1356

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),1356,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1357

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),1357,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1358

ierr = MatSetType(At,((PetscObject)A)->type_name);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

1360

/* get symbolic C=At*B */

1361

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

1362

1363

/* clean up */

1364

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

1365

ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);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

1367

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

1368

1369

}

1370

1371

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

1372

{

1373

PetscErrorCode ierr;

1374

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

1375

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

1376

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

1377

PetscLogDouble flops=0.0;

1378

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

1379

1380

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] = 1380; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1381

if (!c->a) {

1382

ierr = PetscCalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_TRUE,1382,__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),1382,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1383

1384

c->a = ca;

1385

c->free_a = PETSC_TRUE;

1386

} else {

1387

ca = c->a;

1388

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

1389

}

1390

1391

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

1392

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

1393

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

1394

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

1395

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

1396

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

1397

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

1398

nextb = 0;

1399

crow = *aj++;

1400

cjj = cj + ci[crow];

1401

caj = ca + ci[crow];

1402

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

1403

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

1404

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

1405

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

1406

nextb++;

1407

}

1408

}

1409

flops += 2*bnzi;

1410

aa++;

1411

}

1412

}

1413

1414

/* Assemble the final matrix and clean up */

1415

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

1416

ierr = MatAssemblyEnd(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 = PetscLogFlops(flops);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

1419

}

1420

1421

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

1422

{

1423

PetscErrorCode ierr;

1424

1425

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] = 1425; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1426

if (scall == MAT_INITIAL_MATRIX) {

1427

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),1427,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1428

ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);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 = 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),1429,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1430

}

1431

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),1431,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1432

ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);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 = 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),1433,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1434

1435

}

1436

1437

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

1438

{

1439

PetscErrorCode ierr;

1440

1441

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] = 1441; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1442

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

1443

1444

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

1445

1446

}

1447

1448

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

1449

{

1450

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

1451

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

1452

PetscErrorCode ierr;

1453

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

1454

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

1455

const PetscInt *aj;

1456

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

1457

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

1458

1459

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] = 1459; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1460

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

1461

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),1461,__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);

1462

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),1462,__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);

1463

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),1463,__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);

1464

b = bd->v;

1465

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

1466

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

1467

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

1468

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

1469

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

1470

r1 = r2 = r3 = r4 = 0.0;

1471

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

1472

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

1473

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

1474

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

1475

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

1476

r1 += aatmp*b1[ajtmp];

1477

r2 += aatmp*b2[ajtmp];

1478

r3 += aatmp*b3[ajtmp];

1479

r4 += aatmp*b4[ajtmp];

1480

}

1481

c1[i] = r1;

1482

c2[i] = r2;

1483

c3[i] = r3;

1484

c4[i] = r4;

1485

}

1486

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

1487

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

1488

}

1489

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

1490

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

1491

r1 = 0.0;

1492

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

1493

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

1494

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

1495

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

1496

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

1497

}

1498

c1[i] = r1;

1499

}

1500

b1 += bm;

1501

c1 += am;

1502

}

1503

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

1504

ierr = MatDenseRestoreArray(C,&c);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 = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);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 = MatAssemblyEnd(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

1508

}

1509

1510

1511

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

1512

1513

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

1514

{

1515

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

1516

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

1517

PetscErrorCode ierr;

1518

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

1519

MatScalar *aa;

1520

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

1521

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

1522

1523

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] = 1523; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1524

1525

b = bd->v;

1526

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

1527

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

1528

1529

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

1530

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

1531

colam = col*am;

1532

arm = a->compressedrow.nrows;

1533

ii = a->compressedrow.i;

1534

ridx = a->compressedrow.rindex;

1535

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

1536

r1 = r2 = r3 = r4 = 0.0;

1537

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

1538

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

1539

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

1540

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

1541

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

1542

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

1543

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

1544

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

1545

}

1546

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

1547

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

1548

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

1549

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

1550

}

1551

b1 += bm4;

1552

b2 += bm4;

1553

b3 += bm4;

1554

b4 += bm4;

1555

}

1556

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

1557

colam = col*am;

1558

arm = a->compressedrow.nrows;

1559

ii = a->compressedrow.i;

1560

ridx = a->compressedrow.rindex;

1561

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

1562

r1 = 0.0;

1563

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

1564

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

1565

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

1566

1567

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

1568

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

1569

}

1570

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

1571

}

1572

b1 += bm;

1573

}

1574

} else {

1575

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

1576

colam = col*am;

1577

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

1578

r1 = r2 = r3 = r4 = 0.0;

1579

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

1580

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

1581

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

1582

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

1583

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

1584

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

1585

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

1586

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

1587

}

1588

c[colam + i] += r1;

1589

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

1590

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

1591

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

1592

}

1593

b1 += bm4;

1594

b2 += bm4;

1595

b3 += bm4;

1596

b4 += bm4;

1597

}

1598

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

1599

colam = col*am;

1600

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

1601

r1 = 0.0;

1602

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

1603

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

1604

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

1605

1606

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

1607

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

1608

}

1609

c[colam + i] += r1;

1610

}

1611

b1 += bm;

1612

}

1613

}

1614

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

1615

ierr = MatDenseRestoreArray(C,&c);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

1617

}

1618

1619

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

1620

{

1621

PetscErrorCode ierr;

1622

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

1623

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

1624

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

1625

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

1626

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

1627

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

1628

1629

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] = 1629; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1630

btval_den=btdense->v;

1631

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),1631,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1632

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

1633

ncolumns = coloring->ncolumns[k];

1634

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

1635

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

1636

btcol = bj + bi[col];

1637

btval = ba + bi[col];

1638

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

1639

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

1640

brow = btcol[j];

1641

btval_den[brow] = btval[j];

1642

}

1643

}

1644

btval_den += m;

1645

}

1646

1647

}

1648

1649

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

1650

{

1651

PetscErrorCode ierr;

1652

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

1653

const PetscScalar *ca_den,*ca_den_ptr;

1654

PetscScalar *ca=csp->a;

1655

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

1656

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

1657

PetscInt nrows,*row,*idx;

1658

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

1659

1660

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] = 1660; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1661

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

1662

1663

if (brows > 0) {

1664

PetscInt *lstart,row_end,row_start;

1665

lstart = matcoloring->lstart;

1666

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

1667

1668

row_end = brows;

1669

if (row_end > m) row_end = m;

1670

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

1671

ca_den_ptr = ca_den;

1672

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

1673

nrows = matcoloring->nrows[k];

1674

row = rows + colorforrow[k];

1675

idx = den2sp + colorforrow[k];

1676

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

1677

if (row[l] >= row_end) {

1678

lstart[k] = l;

1679

break;

1680

} else {

1681

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

1682

}

1683

}

1684

ca_den_ptr += m;

1685

}

1686

row_end += brows;

1687

if (row_end > m) row_end = m;

1688

}

1689

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

1690

ca_den_ptr = ca_den;

1691

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

1692

nrows = matcoloring->nrows[k];

1693

row = rows + colorforrow[k];

1694

idx = den2sp + colorforrow[k];

1695

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

1696

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

1697

}

1698

ca_den_ptr += m;

1699

}

1700

}

1701

1702

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

1703

#if defined(PETSC_USE_INFO1)

1704

if (matcoloring->brows > 0) {

1705

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),1705,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1706

} else {

1707

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),1707,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1708

}

1709

#endif

1710

1711

}

1712

1713

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

1714

{

1715

PetscErrorCode ierr;

1716

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

1717

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

1718

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

1719

IS *isa;

1720

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

1721

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

1722

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

1723

PetscBool flg;

1724

1725

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] = 1725; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1726

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

1727

1728

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

1729

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

1730

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

1731

c->N = Nbs;

1732

c->m = c->M;

1733

c->rstart = 0;

1734

c->brows = 100;

1735

1736

c->ncolors = nis;

1737

ierr = PetscMalloc3(nis,&c->ncolumns,nis,&c->nrows,nis+1,&colorforrow)PetscMallocA(3,PETSC_FALSE,1737,__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),1737,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1738

ierr = PetscMalloc1(csp->nz+1,&rows)PetscMallocA(1,PETSC_FALSE,1738,__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),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,&den2sp)PetscMallocA(1,PETSC_FALSE,1739,__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),1739,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1740

1741

brows = c->brows;

1742

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),1742,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1743

if (flg) c->brows = brows;

1744

if (brows > 0) {

1745

ierr = PetscMalloc1(nis+1,&c->lstart)PetscMallocA(1,PETSC_FALSE,1745,__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),1745,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1746

}

1747

1748

colorforrow[0] = 0;

1749

rows_i = rows;

1750

den2sp_i = den2sp;

1751

1752

ierr = PetscMalloc1(nis+1,&colorforcol)PetscMallocA(1,PETSC_FALSE,1752,__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),1752,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1753

ierr = PetscMalloc1(Nbs+1,&columns)PetscMallocA(1,PETSC_FALSE,1753,__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),1753,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1754

1755

colorforcol[0] = 0;

1756

columns_i = columns;

1757

1758

/* get column-wise storage of mat */

1759

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),1759,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1760

1761

cm = c->m;

1762

ierr = PetscMalloc1(cm+1,&rowhit)PetscMallocA(1,PETSC_FALSE,1762,__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),1762,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1763

ierr = PetscMalloc1(cm+1,&idxhit)PetscMallocA(1,PETSC_FALSE,1763,__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),1763,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1764

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

1765

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

1766

ierr = ISGetIndices(isa[i],&is);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

1768

c->ncolumns[i] = n;

1769

if (n) {

1770

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),1770,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1771

}

1772

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

1773

columns_i += n;

1774

1775

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

1776

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

1777

1778

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

1779

col = is[j];

1780

row_idx = cj + ci[col];

1781

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

1782

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

1783

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

1784

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

1785

}

1786

}

1787

/* count the number of hits */

1788

nrows = 0;

1789

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

1790

if (rowhit[j]) nrows++;

1791

}

1792

c->nrows[i] = nrows;

1793

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

1794

1795

nrows = 0;

1796

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

1797

if (rowhit[j]) {

1798

rows_i[nrows] = j;

1799

den2sp_i[nrows] = idxhit[j];

1800

nrows++;

1801

}

1802

}

1803

den2sp_i += nrows;

1804

1805

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

1806

rows_i += nrows;

1807

}

1808

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),1808,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1809

ierr = PetscFree(rowhit)((*PetscTrFree)((void*)(rowhit),1809,__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),1809,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1810

ierr = ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&isa);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

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),1811,__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]);

1812

1813

c->colorforrow = colorforrow;

1814

c->rows = rows;

1815

c->den2sp = den2sp;

1816

c->colorforcol = colorforcol;

1817

c->columns = columns;

1818

1819

ierr = PetscFree(idxhit)((*PetscTrFree)((void*)(idxhit),1819,__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),1819,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1820

1821

}

1822

1823

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

1824

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

1825

{

1826

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] = 1826; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1827

/* Empty function */

1828

1829

}

1830

1831

/* This algorithm combines the symbolic and numeric phase of matrix-matrix multiplication. */

1832

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

1833

{

1834

PetscErrorCode ierr;

1835

PetscLogDouble flops=0.0;

1836

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

1837

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

1838

PetscInt *ci,*cj,*cj_i;

1839

PetscScalar *ca,*ca_i;

1840

PetscInt b_maxmemrow,c_maxmem,a_col;

1841

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

1842

PetscInt i,k,ndouble=0;

1843

PetscReal afill;

1844

PetscScalar *c_row_val_dense;

1845

PetscBool *c_row_idx_flags;

1846

PetscInt *aj_i=a->j;

1847

PetscScalar *aa_i=a->a;

1848

1849

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] = 1849; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);

1850

1851

/* Step 1: Determine upper bounds on memory for C and allocate memory */

1852

/* This should be enough for almost all matrices. If still more memory is needed, it is reallocated later. */

1853

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

1854

b_maxmemrow = PetscMin(bi[bm],bn)(((bi[bm])<(bn)) ? (bi[bm]) : (bn));

1855

ierr = PetscMalloc1(am+1,&ci)PetscMallocA(1,PETSC_FALSE,1855,__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),1855,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1856

ierr = PetscCalloc1(bn,&c_row_val_dense)PetscMallocA(1,PETSC_TRUE,1856,__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),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_idx_flags)PetscMallocA(1,PETSC_TRUE,1857,__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),1857,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1858

ierr = PetscMalloc1(c_maxmem,&cj)PetscMallocA(1,PETSC_FALSE,1858,__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),1858,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1859

ierr = PetscMalloc1(c_maxmem,&ca)PetscMallocA(1,PETSC_FALSE,1859,__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),1859,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1860

ca_i = ca;

1861

cj_i = cj;

1862

ci[0] = 0;

1863

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

1864

/* Step 2: Initialize the dense row vector for C */

1865

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

1866

PetscInt cnzi = 0;

1867

PetscInt *bj_i;

1868

PetscScalar *ba_i;

1869

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

1870

Usually, there is enough memory in the first place, so this is not executed. */

1871

while (ci[i] + b_maxmemrow > c_maxmem) {

1872

c_maxmem *= 2;

1873

ndouble++;

1874

ierr = PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj)((*PetscTrRealloc)((sizeof(PetscInt)*c_maxmem),1874,__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),1874,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1875

ierr = PetscRealloc(sizeof(PetscScalar)*c_maxmem,&ca)((*PetscTrRealloc)((sizeof(PetscScalar)*c_maxmem),1875,__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),1875,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1876

}

1877

1878

/* Step 3: Do the numerical calculations */

1879

for (a_col=0; a_col<anzi; a_col++) { /* iterate over all non zero values in a row of A */

1880

PetscInt a_col_index = aj_i[a_col];

1881

const PetscInt bnzi = bi[a_col_index+1] - bi[a_col_index];

1882

flops += 2*bnzi;

1883

bj_i = b->j + bi[a_col_index]; /* points to the current row in bj */

1884

ba_i = b->a + bi[a_col_index]; /* points to the current row in ba */

1885

for (k=0; k<bnzi; ++k) { /* iterate over all non zeros of this row in B */

1886

if (c_row_idx_flags[bj_i[k]] == PETSC_FALSE) {

1887

cj_i[cnzi++] = bj_i[k];

1888

c_row_idx_flags[bj_i[k]] = PETSC_TRUE;

1889

}

1890

c_row_val_dense[bj_i[k]] += aa_i[a_col] * ba_i[k];

1891

}

1892

}

1893

1894

/* Sort array */

1895

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

1896

/* Step 4 */

1897

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

1898

ca_i[k] = c_row_val_dense[cj_i[k]];

1899

c_row_val_dense[cj_i[k]] = 0.;

1900

c_row_idx_flags[cj_i[k]] = PETSC_FALSE;

1901

}

1902

/* terminate current row */

1903

aa_i += anzi;

1904

aj_i += anzi;

1905

ca_i += cnzi;

1906

cj_i += cnzi;

1907

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

1908

flops += cnzi;

1909

}

1910

1911

/* Step 5 */

1912

/* Create the new matrix */

1913

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),1913,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1914

ierr = MatSetBlockSizesFromMats(*C,A,B);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 = MatSetType(*C,((PetscObject)A)->type_name);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

1917

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

1918

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

1919

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

1920

c->a = ca;

1921

c->free_a = PETSC_TRUE;

1922

c->free_ij = PETSC_TRUE;

1923

c->nonew = 0;

1924

1925

/* set MatInfo */

1926

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

1927

if (afill < 1.0) afill = 1.0;

1928

c->maxnz = ci[am];

1929

c->nz = ci[am];

1930

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

1931

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

1932

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

1933

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

1934

1935

ierr = PetscFree(c_row_val_dense)((*PetscTrFree)((void*)(c_row_val_dense),1935,__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),1935,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1936

ierr = PetscFree(c_row_idx_flags)((*PetscTrFree)((void*)(c_row_idx_flags),1936,__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),1936,__func__,"/sandbox/petsc/petsc.next/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);

1937

1938

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

1939

ierr = MatAssemblyEnd(*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 = PetscLogFlops(flops);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

1942

}