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

1
2/*
3 Defines matrix-matrix product routines for pairs of SeqAIJ matrices
4 C = A * B
5*/
6
7#include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
8#include <../src/mat/utils/freespace.h>
9#include <petscbt.h>
10#include <petsc/private/isimpl.h>
11#include <../src/mat/impls/dense/seq/dense.h>
12
13
14PETSC_INTERNextern __attribute__((visibility ("hidden"))) PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
15{
16 PetscErrorCode ierr;
17
18 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)
;
19 if (scall == MAT_INITIAL_MATRIX) {
20 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)
;
21 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)
;
22 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)
;
23 }
24 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)
;
25 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)
;
26 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)
;
27 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)
;
28}
29
30PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
31{
32 PetscErrorCode ierr;
33
34 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)
;
35 if (C->ops->matmultnumeric) {
36 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)
;
37 } else {
38 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)
;
39 }
40 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)
;
41}
42
43PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
44{
45 PetscErrorCode ierr;
46#if !defined(PETSC_HAVE_HYPRE)
47 const char *algTypes[8] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge"};
48 PetscInt nalg = 8;
49#else
50 const char *algTypes[9] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge","hypre"};
51 PetscInt nalg = 9;
52#endif
53 PetscInt alg = 0; /* set default algorithm */
54
55 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)
;
56 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)
;
57 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)
;
58 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)
;
59 switch (alg) {
60 case 1:
61 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)
;
62 break;
63 case 2:
64 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)
;
65 break;
66 case 3:
67 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)
;
68 break;
69 case 4:
70 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)
;
71 break;
72 case 5:
73 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)
;
74 break;
75 case 6:
76 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)
;
77 break;
78 case 7:
79 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)
;
80 break;
81#if defined(PETSC_HAVE_HYPRE)
82 case 8:
83 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)
;
84 break;
85#endif
86 default:
87 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)
;
88 break;
89 }
90 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)
;
91}
92
93PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(Mat A,Mat B,PetscReal fill,Mat *C)
94{
95 PetscErrorCode ierr;
96 Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
97 PetscInt *ai=a->i,*bi=b->i,*ci,*cj;
98 PetscInt am =A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
99 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 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)
;
194}
195
196PetscErrorCode 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 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)
;
257}
258
259PetscErrorCode 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 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)
;
310}
311
312PetscErrorCode 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 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); } } } }
;
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 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)
;
415}
416
417
418PetscErrorCode 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 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); } } } }
;
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 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)
;
521}
522
523PetscErrorCode 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 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)
;
628}
629
630PetscErrorCode 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;
1
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++) {
2
Assuming 'i' is < 'am'
3
Loop condition is true. Entering loop body
11
Assuming 'i' is < 'am'
12
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 */
4
Assuming pointer value is null
664 ci[i+1] = ci[i];
665 /* Populate the min heap */
666 for (j=0; j<anzi; j++) {
5
Assuming 'j' is >= 'anzi'
6
Loop condition is false. Execution continues on line 678
13
Assuming 'j' is >= 'anzi'
14
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) {
7
Assuming 'j' is < 0
8
Loop condition is false. Execution continues on line 701
15
Assuming 'j' is >= 0
16
Loop condition is true. Entering loop body
680 if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
17
Assuming the condition is false
18
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;
19
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 */
9
Assuming 'fptr' is null
10
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 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)
;
749}
750
751
752PetscErrorCode 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 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)
;
1014}
1015
1016/* concatenate unique entries and then sort */
1017PetscErrorCode 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 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)
;
1102}
1103
1104/* This routine is not used. Should be removed! */
1105PetscErrorCode 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 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)
;
1119}
1120
1121PetscErrorCode 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 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)
;
1134}
1135
1136PetscErrorCode 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 PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0)
;
1214}
1215
1216PetscErrorCode 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 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)
;
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 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)
;
1285}
1286
1287PetscErrorCode 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 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)
;
1300}
1301
1302PetscErrorCode 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 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)
;
1345}
1346
1347PetscErrorCode 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 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)
;
1369}
1370
1371PetscErrorCode 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 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)
;
1419}
1420
1421PETSC_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 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)
;
1435}
1436
1437PetscErrorCode 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 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)
;
1446}
1447
1448PetscErrorCode 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 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)
;
1508}
1509
1510/*
1511 Note very similar to MatMult_SeqAIJ(), should generate both codes from same base
1512*/
1513PetscErrorCode 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 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)
;
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 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)
;
1617}
1618
1619PetscErrorCode 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 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)
;
1647}
1648
1649PetscErrorCode 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 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)
;
1711}
1712
1713PetscErrorCode 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 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)
;
1821}
1822
1823/* The combine method has done the symbolic and numeric in the first phase, and so we just return here */
1824PetscErrorCode 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 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)
;
1829}
1830
1831/* This algorithm combines the symbolic and numeric phase of matrix-matrix multiplication. */
1832PetscErrorCode 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 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)
;
1942}