Bug Summary

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

Annotated Source Code

[?] Use j/k keys for keyboard navigation

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
14 PETSC_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-tmp/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-tmp/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-tmp/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-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
23 }
24
25 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultNumeric].active) ? (*PetscLogPLB)((MAT_MatMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),25,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
26 ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),26,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
27 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultNumeric].active) ? (*PetscLogPLE)((MAT_MatMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),27,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
28 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)
;
29 }
30
31 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
32 {
33 PetscErrorCode ierr;
34
35 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 35; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
36 if (C->ops->matmultnumeric) {
37 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),37,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
38 } else {
39 ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(A,B,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),39,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
40 }
41 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)
;
42 }
43
44 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
45 {
46 PetscErrorCode ierr;
47 #if !defined(PETSC_HAVE_HYPRE)
48 const char *algTypes[8] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge"};
49 PetscInt nalg = 8;
50 #else
51 const char *algTypes[9] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge","hypre"};
52 PetscInt nalg = 9;
53 #endif
54 PetscInt alg = 0; /* set default algorithm */
55
56 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 56; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
57 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatMatMult","Mat")0; do { PetscOptionItems PetscOptionsObjectBase; PetscOptionItems
*PetscOptionsObject = &PetscOptionsObjectBase; PetscMemzero
(PetscOptionsObject,sizeof(PetscOptionItems)); for (PetscOptionsObject
->count=(PetscOptionsPublish?-1:1); PetscOptionsObject->
count<2; PetscOptionsObject->count++) { PetscErrorCode _5_ierr
= PetscOptionsBegin_Private(PetscOptionsObject,PetscObjectComm
((PetscObject)A),((PetscObject)A)->prefix,"MatMatMult","Mat"
);do {if (__builtin_expect(!!(_5_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),57,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),57,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
58 ierr = PetscOptionsEList("-matmatmult_via","Algorithmic approach","MatMatMult",algTypes,nalg,algTypes[0],&alg,NULL)PetscOptionsEList_Private(PetscOptionsObject,"-matmatmult_via"
,"Algorithmic approach","MatMatMult",algTypes,nalg,algTypes[0
],&alg,((void*)0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),58,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
59 ierr = PetscOptionsEnd()_5_ierr = PetscOptionsEnd_Private(PetscOptionsObject);do {if (
__builtin_expect(!!(_5_ierr),0)) return PetscError(((MPI_Comm
)0x44000001),59,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);}} while (0)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),59,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
60 switch (alg) {
61 case 1:
62 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),62,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
63 break;
64 case 2:
65 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),65,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
66 break;
67 case 3:
68 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),68,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
69 break;
70 case 4:
71 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),71,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
72 break;
73 case 5:
74 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),74,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
75 break;
76 case 6:
77 ierr = MatMatMult_SeqAIJ_SeqAIJ_Combined(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),77,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
78 break;
79 case 7:
80 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),80,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
81 break;
82 #if defined(PETSC_HAVE_HYPRE)
83 case 8:
84 ierr = MatMatMultSymbolic_AIJ_AIJ_wHYPRE(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),84,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
85 break;
86 #endif
87 default:
88 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),88,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
89 break;
90 }
91 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)
;
92 }
93
94 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(Mat A,Mat B,PetscReal fill,Mat *C)
95 {
96 PetscErrorCode ierr;
97 Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
98 PetscInt *ai=a->i,*bi=b->i,*ci,*cj;
99 PetscInt am =A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
100 PetscReal afill;
101 PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
102 PetscTable ta;
103 PetscBT lnkbt;
104 PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);
105
106 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 106; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
107 /* Get ci and cj */
108 /*---------------*/
109 /* Allocate ci array, arrays for fill computation and */
110 /* free space for accumulating nonzero column info */
111 ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,111,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),111,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
112 ci[0] = 0;
113
114 /* create and initialize a linked list */
115 ierr = PetscTableCreate(bn,bn,&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),115,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
116 MatRowMergeMax_SeqAIJ(b,bm,ta){ PetscInt _j,_row,_nz,*_col; if (b) { for (_row=0; _row<bm
; _row++) { _nz = b->i[_row+1] - b->i[_row]; for (_j=0;
_j<_nz; _j++) { _col = _j + b->j + b->i[_row]; PetscTableAdd
(ta,*_col+1,1,INSERT_VALUES); } } } }
;
117 ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),117,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
118 ierr = PetscTableDestroy(&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),118,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
119
120 ierr = PetscLLCondensedCreate(Crmax,bn,&lnk,&lnkbt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),120,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
121
122 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
123 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),123,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
124
125 current_space = free_space;
126
127 /* Determine ci and cj */
128 for (i=0; i<am; i++) {
129 anzi = ai[i+1] - ai[i];
130 aj = a->j + ai[i];
131 for (j=0; j<anzi; j++) {
132 brow = aj[j];
133 bnzj = bi[brow+1] - bi[brow];
134 bj = b->j + bi[brow];
135 /* add non-zero cols of B into the sorted linked list lnk */
136 ierr = PetscLLCondensedAddSorted(bnzj,bj,lnk,lnkbt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),136,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
137 }
138 cnzi = lnk[0];
139
140 /* If free space is not available, make more free space */
141 /* Double the amount of total space in the list */
142 if (current_space->local_remaining<cnzi) {
143 ierr = PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),143,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
144 ndouble++;
145 }
146
147 /* Copy data into free space, then initialize lnk */
148 ierr = PetscLLCondensedClean(bn,cnzi,current_space->array,lnk,lnkbt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),148,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
149
150 current_space->array += cnzi;
151 current_space->local_used += cnzi;
152 current_space->local_remaining -= cnzi;
153
154 ci[i+1] = ci[i] + cnzi;
155 }
156
157 /* Column indices are in the list of free space */
158 /* Allocate space for cj, initialize cj, and */
159 /* destroy list of free space and other temporary array(s) */
160 ierr = PetscMalloc1(ci[am]+1,&cj)PetscMallocA(1,PETSC_FALSE,160,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&cj)),(&cj))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),160,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
161 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),161,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
162 ierr = PetscLLCondensedDestroy(lnk,lnkbt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),162,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
163
164 /* put together the new symbolic matrix */
165 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),165,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
166 ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),166,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
167 ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),167,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
168
169 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
170 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
171 c = (Mat_SeqAIJ*)((*C)->data);
172 c->free_a = PETSC_FALSE;
173 c->free_ij = PETSC_TRUE;
174 c->nonew = 0;
175 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted; /* fast, needs non-scalable O(bn) array 'abdense' */
176
177 /* set MatInfo */
178 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
179 if (afill < 1.0) afill = 1.0;
180 c->maxnz = ci[am];
181 c->nz = ci[am];
182 (*C)->info.mallocs = ndouble;
183 (*C)->info.fill_ratio_given = fill;
184 (*C)->info.fill_ratio_needed = afill;
185
186#if defined(PETSC_USE_INFO1)
187 if (ci[am]) {
188 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),188,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
189 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),189,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
190 } else {
191 ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),191,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
192 }
193#endif
194 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)
;
195}
196
197PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,Mat C)
198{
199 PetscErrorCode ierr;
200 PetscLogDouble flops=0.0;
201 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
202 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
203 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
204 PetscInt *ai =a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;
205 PetscInt am =A->rmap->n,cm=C->rmap->n;
206 PetscInt i,j,k,anzi,bnzi,cnzi,brow;
207 PetscScalar *aa=a->a,*ba=b->a,*baj,*ca,valtmp;
208 PetscScalar *ab_dense;
209
210 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 210; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
211 if (!c->a) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
212 ierr = PetscMalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_FALSE,212,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[cm]+1)*sizeof(**(&ca)),(&ca))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),212,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
213 c->a = ca;
214 c->free_a = PETSC_TRUE;
215 } else {
216 ca = c->a;
217 }
218 if (!c->matmult_abdense) {
219 ierr = PetscCalloc1(B->cmap->N,&ab_dense)PetscMallocA(1,PETSC_TRUE,219,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(B->cmap->N)*sizeof(**(&ab_dense)),(&ab_dense
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),219,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
220 c->matmult_abdense = ab_dense;
221 } else {
222 ab_dense = c->matmult_abdense;
223 }
224
225 /* clean old values in C */
226 ierr = PetscArrayzero(ca,ci[cm])PetscMemzero(ca,(ci[cm])*sizeof(*(ca)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),226,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
227 /* Traverse A row-wise. */
228 /* Build the ith row in C by summing over nonzero columns in A, */
229 /* the rows of B corresponding to nonzeros of A. */
230 for (i=0; i<am; i++) {
231 anzi = ai[i+1] - ai[i];
232 for (j=0; j<anzi; j++) {
233 brow = aj[j];
234 bnzi = bi[brow+1] - bi[brow];
235 bjj = bj + bi[brow];
236 baj = ba + bi[brow];
237 /* perform dense axpy */
238 valtmp = aa[j];
239 for (k=0; k<bnzi; k++) {
240 ab_dense[bjj[k]] += valtmp*baj[k];
241 }
242 flops += 2*bnzi;
243 }
244 aj += anzi; aa += anzi;
245
246 cnzi = ci[i+1] - ci[i];
247 for (k=0; k<cnzi; k++) {
248 ca[k] += ab_dense[cj[k]];
249 ab_dense[cj[k]] = 0.0; /* zero ab_dense */
250 }
251 flops += cnzi;
252 cj += cnzi; ca += cnzi;
253 }
254 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),254,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
255 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),255,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
256 ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),256,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
257 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)
;
258}
259
260PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,Mat C)
261{
262 PetscErrorCode ierr;
263 PetscLogDouble flops=0.0;
264 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
265 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
266 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
267 PetscInt *ai = a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;
268 PetscInt am = A->rmap->N,cm=C->rmap->N;
269 PetscInt i,j,k,anzi,bnzi,cnzi,brow;
270 PetscScalar *aa=a->a,*ba=b->a,*baj,*ca=c->a,valtmp;
271 PetscInt nextb;
272
273 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 273; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
274 if (!ca) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
275 ierr = PetscMalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_FALSE,275,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[cm]+1)*sizeof(**(&ca)),(&ca))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),275,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
276 c->a = ca;
277 c->free_a = PETSC_TRUE;
278 }
279
280 /* clean old values in C */
281 ierr = PetscArrayzero(ca,ci[cm])PetscMemzero(ca,(ci[cm])*sizeof(*(ca)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),281,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
282 /* Traverse A row-wise. */
283 /* Build the ith row in C by summing over nonzero columns in A, */
284 /* the rows of B corresponding to nonzeros of A. */
285 for (i=0; i<am; i++) {
286 anzi = ai[i+1] - ai[i];
287 cnzi = ci[i+1] - ci[i];
288 for (j=0; j<anzi; j++) {
289 brow = aj[j];
290 bnzi = bi[brow+1] - bi[brow];
291 bjj = bj + bi[brow];
292 baj = ba + bi[brow];
293 /* perform sparse axpy */
294 valtmp = aa[j];
295 nextb = 0;
296 for (k=0; nextb<bnzi; k++) {
297 if (cj[k] == bjj[nextb]) { /* ccol == bcol */
298 ca[k] += valtmp*baj[nextb++];
299 }
300 }
301 flops += 2*bnzi;
302 }
303 aj += anzi; aa += anzi;
304 cj += cnzi; ca += cnzi;
305 }
306
307 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),307,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
308 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),308,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
309 ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),309,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
310 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)
;
311}
312
313PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(Mat A,Mat B,PetscReal fill,Mat *C)
314{
315 PetscErrorCode ierr;
316 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
317 PetscInt *ai = a->i,*bi=b->i,*ci,*cj;
318 PetscInt am = A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
319 MatScalar *ca;
320 PetscReal afill;
321 PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
322 PetscTable ta;
323 PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);
324
325 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 325; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
326 /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_fast() */
327 /*-----------------------------------------------------------------------------------------*/
328 /* Allocate arrays for fill computation and free space for accumulating nonzero column */
329 ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,329,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),329,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
330 ci[0] = 0;
331
332 /* create and initialize a linked list */
333 ierr = PetscTableCreate(bn,bn,&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),333,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
334 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); } } } }
;
335 ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),335,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
336 ierr = PetscTableDestroy(&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),336,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
337
338 ierr = PetscLLCondensedCreate_fast(Crmax,&lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),338,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
339
340 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
341 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),341,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
342 current_space = free_space;
343
344 /* Determine ci and cj */
345 for (i=0; i<am; i++) {
346 anzi = ai[i+1] - ai[i];
347 aj = a->j + ai[i];
348 for (j=0; j<anzi; j++) {
349 brow = aj[j];
350 bnzj = bi[brow+1] - bi[brow];
351 bj = b->j + bi[brow];
352 /* add non-zero cols of B into the sorted linked list lnk */
353 ierr = PetscLLCondensedAddSorted_fast(bnzj,bj,lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),353,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
354 }
355 cnzi = lnk[1];
356
357 /* If free space is not available, make more free space */
358 /* Double the amount of total space in the list */
359 if (current_space->local_remaining<cnzi) {
360 ierr = PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),360,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
361 ndouble++;
362 }
363
364 /* Copy data into free space, then initialize lnk */
365 ierr = PetscLLCondensedClean_fast(cnzi,current_space->array,lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),365,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
366
367 current_space->array += cnzi;
368 current_space->local_used += cnzi;
369 current_space->local_remaining -= cnzi;
370
371 ci[i+1] = ci[i] + cnzi;
372 }
373
374 /* Column indices are in the list of free space */
375 /* Allocate space for cj, initialize cj, and */
376 /* destroy list of free space and other temporary array(s) */
377 ierr = PetscMalloc1(ci[am]+1,&cj)PetscMallocA(1,PETSC_FALSE,377,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&cj)),(&cj))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),377,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
378 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),378,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
379 ierr = PetscLLCondensedDestroy_fast(lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),379,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
380
381 /* Allocate space for ca */
382 ierr = PetscCalloc1(ci[am]+1,&ca)PetscMallocA(1,PETSC_TRUE,382,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&ca)),(&ca))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),382,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
383
384 /* put together the new symbolic matrix */
385 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,ca,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),385,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
386 ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),386,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
387 ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),387,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
388
389 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
390 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
391 c = (Mat_SeqAIJ*)((*C)->data);
392 c->free_a = PETSC_TRUE;
393 c->free_ij = PETSC_TRUE;
394 c->nonew = 0;
395
396 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */
397
398 /* set MatInfo */
399 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
400 if (afill < 1.0) afill = 1.0;
401 c->maxnz = ci[am];
402 c->nz = ci[am];
403 (*C)->info.mallocs = ndouble;
404 (*C)->info.fill_ratio_given = fill;
405 (*C)->info.fill_ratio_needed = afill;
406
407#if defined(PETSC_USE_INFO1)
408 if (ci[am]) {
409 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),409,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
410 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),410,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
411 } else {
412 ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),412,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
413 }
414#endif
415 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)
;
416}
417
418
419PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,PetscReal fill,Mat *C)
420{
421 PetscErrorCode ierr;
422 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
423 PetscInt *ai = a->i,*bi=b->i,*ci,*cj;
424 PetscInt am = A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
425 MatScalar *ca;
426 PetscReal afill;
427 PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
428 PetscTable ta;
429 PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);
430
431 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 431; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
432 /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_Scalalbe() */
433 /*---------------------------------------------------------------------------------------------*/
434 /* Allocate arrays for fill computation and free space for accumulating nonzero column */
435 ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,435,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),435,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
436 ci[0] = 0;
437
438 /* create and initialize a linked list */
439 ierr = PetscTableCreate(bn,bn,&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),439,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
440 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); } } } }
;
441 ierr = PetscTableGetCount(ta,&Crmax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),441,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
442 ierr = PetscTableDestroy(&ta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),442,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
443 ierr = PetscLLCondensedCreate_Scalable(Crmax,&lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),443,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
444
445 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
446 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),446,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
447 current_space = free_space;
448
449 /* Determine ci and cj */
450 for (i=0; i<am; i++) {
451 anzi = ai[i+1] - ai[i];
452 aj = a->j + ai[i];
453 for (j=0; j<anzi; j++) {
454 brow = aj[j];
455 bnzj = bi[brow+1] - bi[brow];
456 bj = b->j + bi[brow];
457 /* add non-zero cols of B into the sorted linked list lnk */
458 ierr = PetscLLCondensedAddSorted_Scalable(bnzj,bj,lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),458,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
459 }
460 cnzi = lnk[0];
461
462 /* If free space is not available, make more free space */
463 /* Double the amount of total space in the list */
464 if (current_space->local_remaining<cnzi) {
465 ierr = PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),465,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
466 ndouble++;
467 }
468
469 /* Copy data into free space, then initialize lnk */
470 ierr = PetscLLCondensedClean_Scalable(cnzi,current_space->array,lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),470,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
471
472 current_space->array += cnzi;
473 current_space->local_used += cnzi;
474 current_space->local_remaining -= cnzi;
475
476 ci[i+1] = ci[i] + cnzi;
477 }
478
479 /* Column indices are in the list of free space */
480 /* Allocate space for cj, initialize cj, and */
481 /* destroy list of free space and other temporary array(s) */
482 ierr = PetscMalloc1(ci[am]+1,&cj)PetscMallocA(1,PETSC_FALSE,482,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&cj)),(&cj))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),482,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
483 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),483,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
484 ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),484,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
485
486 /* Allocate space for ca */
487 /*-----------------------*/
488 ierr = PetscCalloc1(ci[am]+1,&ca)PetscMallocA(1,PETSC_TRUE,488,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am]+1)*sizeof(**(&ca)),(&ca))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),488,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
489
490 /* put together the new symbolic matrix */
491 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,ca,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),491,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
492 ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),492,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
493 ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),493,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
494
495 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
496 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
497 c = (Mat_SeqAIJ*)((*C)->data);
498 c->free_a = PETSC_TRUE;
499 c->free_ij = PETSC_TRUE;
500 c->nonew = 0;
501
502 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */
503
504 /* set MatInfo */
505 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
506 if (afill < 1.0) afill = 1.0;
507 c->maxnz = ci[am];
508 c->nz = ci[am];
509 (*C)->info.mallocs = ndouble;
510 (*C)->info.fill_ratio_given = fill;
511 (*C)->info.fill_ratio_needed = afill;
512
513#if defined(PETSC_USE_INFO1)
514 if (ci[am]) {
515 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),515,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
516 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),516,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
517 } else {
518 ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),518,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
519 }
520#endif
521 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)
;
522}
523
524PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(Mat A,Mat B,PetscReal fill,Mat *C)
525{
526 PetscErrorCode ierr;
527 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
528 const PetscInt *ai=a->i,*bi=b->i,*aj=a->j,*bj=b->j;
529 PetscInt *ci,*cj,*bb;
530 PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
531 PetscReal afill;
532 PetscInt i,j,col,ndouble = 0;
533 PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);
534 PetscHeap h;
535
536 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 536; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
537 /* Get ci and cj - by merging sorted rows using a heap */
538 /*---------------------------------------------------------------------------------------------*/
539 /* Allocate arrays for fill computation and free space for accumulating nonzero column */
540 ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,540,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),540,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
541 ci[0] = 0;
542
543 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
544 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),544,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
545 current_space = free_space;
546
547 ierr = PetscHeapCreate(a->rmax,&h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),547,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
548 ierr = PetscMalloc1(a->rmax,&bb)PetscMallocA(1,PETSC_FALSE,548,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a->rmax)*sizeof(**(&bb)),(&bb))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),548,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
549
550 /* Determine ci and cj */
551 for (i=0; i<am; i++) {
552 const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
553 const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
554 ci[i+1] = ci[i];
555 /* Populate the min heap */
556 for (j=0; j<anzi; j++) {
557 bb[j] = bi[acol[j]]; /* bb points at the start of the row */
558 if (bb[j] < bi[acol[j]+1]) { /* Add if row is nonempty */
559 ierr = PetscHeapAdd(h,j,bj[bb[j]++]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),559,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
560 }
561 }
562 /* Pick off the min element, adding it to free space */
563 ierr = PetscHeapPop(h,&j,&col);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),563,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
564 while (j >= 0) {
565 if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
566 ierr = PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20)(((PetscIntMultTruncate(2,current_space->total_array_size)
)<(16 << 20)) ? (PetscIntMultTruncate(2,current_space
->total_array_size)) : (16 << 20))
,&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),566,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
567 ndouble++;
568 }
569 *(current_space->array++) = col;
570 current_space->local_used++;
571 current_space->local_remaining--;
572 ci[i+1]++;
573
574 /* stash if anything else remains in this row of B */
575 if (bb[j] < bi[acol[j]+1]) {ierr = PetscHeapStash(h,j,bj[bb[j]++]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),575,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;}
576 while (1) { /* pop and stash any other rows of B that also had an entry in this column */
577 PetscInt j2,col2;
578 ierr = PetscHeapPeek(h,&j2,&col2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),578,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
579 if (col2 != col) break;
580 ierr = PetscHeapPop(h,&j2,&col2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),580,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
581 if (bb[j2] < bi[acol[j2]+1]) {ierr = PetscHeapStash(h,j2,bj[bb[j2]++]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),581,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;}
582 }
583 /* Put any stashed elements back into the min heap */
584 ierr = PetscHeapUnstash(h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),584,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
585 ierr = PetscHeapPop(h,&j,&col);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),585,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
586 }
587 }
588 ierr = PetscFree(bb)((*PetscTrFree)((void*)(bb),588,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((bb) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),588,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
589 ierr = PetscHeapDestroy(&h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),589,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
590
591 /* Column indices are in the list of free space */
592 /* Allocate space for cj, initialize cj, and */
593 /* destroy list of free space and other temporary array(s) */
594 ierr = PetscMalloc1(ci[am],&cj)PetscMallocA(1,PETSC_FALSE,594,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am])*sizeof(**(&cj)),(&cj))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),594,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
595 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),595,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
596
597 /* put together the new symbolic matrix */
598 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),598,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
599 ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),599,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
600 ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),600,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
601
602 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
603 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
604 c = (Mat_SeqAIJ*)((*C)->data);
605 c->free_a = PETSC_TRUE;
606 c->free_ij = PETSC_TRUE;
607 c->nonew = 0;
608
609 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
610
611 /* set MatInfo */
612 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
613 if (afill < 1.0) afill = 1.0;
614 c->maxnz = ci[am];
615 c->nz = ci[am];
616 (*C)->info.mallocs = ndouble;
617 (*C)->info.fill_ratio_given = fill;
618 (*C)->info.fill_ratio_needed = afill;
619
620#if defined(PETSC_USE_INFO1)
621 if (ci[am]) {
622 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),622,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
623 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),623,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
624 } else {
625 ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),625,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
626 }
627#endif
628 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)
;
629}
630
631PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(Mat A,Mat B,PetscReal fill,Mat *C)
632{
633 PetscErrorCode ierr;
634 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
635 const PetscInt *ai = a->i,*bi=b->i,*aj=a->j,*bj=b->j;
636 PetscInt *ci,*cj,*bb;
637 PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
638 PetscReal afill;
639 PetscInt i,j,col,ndouble = 0;
640 PetscFreeSpaceList free_space=NULL((void*)0),current_space=NULL((void*)0);
641 PetscHeap h;
642 PetscBT bt;
643
644 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 644; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
645 /* Get ci and cj - using a heap for the sorted rows, but use BT so that each index is only added once */
646 /*---------------------------------------------------------------------------------------------*/
647 /* Allocate arrays for fill computation and free space for accumulating nonzero column */
648 ierr = PetscMalloc1(am+2,&ci)PetscMallocA(1,PETSC_FALSE,648,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+2)*sizeof(**(&ci)),(&ci))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),648,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
649 ci[0] = 0;
650
651 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
652 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),652,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
653
654 current_space = free_space;
1
Value assigned to 'current_space'
655
656 ierr = PetscHeapCreate(a->rmax,&h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),656,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
657 ierr = PetscMalloc1(a->rmax,&bb)PetscMallocA(1,PETSC_FALSE,657,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a->rmax)*sizeof(**(&bb)),(&bb))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),657,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
658 ierr = PetscBTCreate(bn,&bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),658,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
659
660 /* Determine ci and cj */
661 for (i=0; i<am; i++) {
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
662 const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
663 const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
664 const PetscInt *fptr = current_space->array; /* Save beginning of the row so we can clear the BT later */
4
Assuming pointer value is null
665 ci[i+1] = ci[i];
666 /* Populate the min heap */
667 for (j=0; j<anzi; j++) {
5
Assuming 'j' is >= 'anzi'
6
Loop condition is false. Execution continues on line 679
13
Assuming 'j' is >= 'anzi'
14
Loop condition is false. Execution continues on line 679
668 PetscInt brow = acol[j];
669 for (bb[j] = bi[brow]; bb[j] < bi[brow+1]; bb[j]++) {
670 PetscInt bcol = bj[bb[j]];
671 if (!PetscBTLookupSet(bt,bcol)) { /* new entry */
672 ierr = PetscHeapAdd(h,j,bcol);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),672,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
673 bb[j]++;
674 break;
675 }
676 }
677 }
678 /* Pick off the min element, adding it to free space */
679 ierr = PetscHeapPop(h,&j,&col);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),679,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
680 while (j >= 0) {
7
Assuming 'j' is < 0
8
Loop condition is false. Execution continues on line 702
15
Assuming 'j' is >= 0
16
Loop condition is true. Entering loop body
681 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
682 fptr = NULL((void*)0); /* need PetscBTMemzero */
683 ierr = PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20)(((PetscIntMultTruncate(2,current_space->total_array_size)
)<(16 << 20)) ? (PetscIntMultTruncate(2,current_space
->total_array_size)) : (16 << 20))
,&current_space);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),683,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
684 ndouble++;
685 }
686 *(current_space->array++) = col;
19
Dereference of null pointer
687 current_space->local_used++;
688 current_space->local_remaining--;
689 ci[i+1]++;
690
691 /* stash if anything else remains in this row of B */
692 for (; bb[j] < bi[acol[j]+1]; bb[j]++) {
693 PetscInt bcol = bj[bb[j]];
694 if (!PetscBTLookupSet(bt,bcol)) { /* new entry */
695 ierr = PetscHeapAdd(h,j,bcol);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),695,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
696 bb[j]++;
697 break;
698 }
699 }
700 ierr = PetscHeapPop(h,&j,&col);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),700,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
701 }
702 if (fptr) { /* Clear the bits for this row */
9
Assuming 'fptr' is null
10
Taking false branch
703 for (; fptr<current_space->array; fptr++) {ierr = PetscBTClear(bt,*fptr);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),703,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;}
704 } else { /* We reallocated so we don't remember (easily) how to clear only the bits we changed */
705 ierr = PetscBTMemzero(bn,bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),705,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
706 }
707 }
708 ierr = PetscFree(bb)((*PetscTrFree)((void*)(bb),708,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((bb) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),708,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
709 ierr = PetscHeapDestroy(&h);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),709,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
710 ierr = PetscBTDestroy(&bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),710,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
711
712 /* Column indices are in the list of free space */
713 /* Allocate space for cj, initialize cj, and */
714 /* destroy list of free space and other temporary array(s) */
715 ierr = PetscMalloc1(ci[am],&cj)PetscMallocA(1,PETSC_FALSE,715,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[am])*sizeof(**(&cj)),(&cj))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),715,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
716 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),716,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
717
718 /* put together the new symbolic matrix */
719 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),719,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
720 ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),720,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
721 ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),721,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
722
723 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
724 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
725 c = (Mat_SeqAIJ*)((*C)->data);
726 c->free_a = PETSC_TRUE;
727 c->free_ij = PETSC_TRUE;
728 c->nonew = 0;
729
730 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
731
732 /* set MatInfo */
733 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
734 if (afill < 1.0) afill = 1.0;
735 c->maxnz = ci[am];
736 c->nz = ci[am];
737 (*C)->info.mallocs = ndouble;
738 (*C)->info.fill_ratio_given = fill;
739 (*C)->info.fill_ratio_needed = afill;
740
741#if defined(PETSC_USE_INFO1)
742 if (ci[am]) {
743 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),743,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
744 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),744,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
745 } else {
746 ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),746,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
747 }
748#endif
749 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)
;
750}
751
752
753PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(Mat A,Mat B,PetscReal fill,Mat *C)
754{
755 PetscErrorCode ierr;
756 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
757 const PetscInt *ai=a->i,*bi=b->i,*aj=a->j,*bj=b->j,*inputi,*inputj,*inputcol,*inputcol_L1;
758 PetscInt *ci,*cj,*outputj,worki_L1[9],worki_L2[9];
759 PetscInt c_maxmem,a_maxrownnz=0,a_rownnz;
760 const PetscInt workcol[8]={0,1,2,3,4,5,6,7};
761 const PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
762 const PetscInt *brow_ptr[8],*brow_end[8];
763 PetscInt window[8];
764 PetscInt window_min,old_window_min,ci_nnz,outputi_nnz=0,L1_nrows,L2_nrows;
765 PetscInt i,k,ndouble=0,L1_rowsleft,rowsleft;
766 PetscReal afill;
767 PetscInt *workj_L1,*workj_L2,*workj_L3;
768 PetscInt L1_nnz,L2_nnz;
769
770 /* Step 1: Get upper bound on memory required for allocation.
771 Because of the way virtual memory works,
772 only the memory pages that are actually needed will be physically allocated. */
773 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 773; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
774 ierr = PetscMalloc1(am+1,&ci)PetscMallocA(1,PETSC_FALSE,774,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+1)*sizeof(**(&ci)),(&ci))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),774,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
775 for (i=0; i<am; i++) {
776 const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
777 const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
778 a_rownnz = 0;
779 for (k=0; k<anzi; ++k) {
780 a_rownnz += bi[acol[k]+1] - bi[acol[k]];
781 if (a_rownnz > bn) {
782 a_rownnz = bn;
783 break;
784 }
785 }
786 a_maxrownnz = PetscMax(a_maxrownnz, a_rownnz)(((a_maxrownnz)<(a_rownnz)) ? (a_rownnz) : (a_maxrownnz));
787 }
788 /* temporary work areas for merging rows */
789 ierr = PetscMalloc1(a_maxrownnz*8,&workj_L1)PetscMallocA(1,PETSC_FALSE,789,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a_maxrownnz*8)*sizeof(**(&workj_L1)),(&workj_L1
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),789,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
790 ierr = PetscMalloc1(a_maxrownnz*8,&workj_L2)PetscMallocA(1,PETSC_FALSE,790,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a_maxrownnz*8)*sizeof(**(&workj_L2)),(&workj_L2
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),790,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
791 ierr = PetscMalloc1(a_maxrownnz,&workj_L3)PetscMallocA(1,PETSC_FALSE,791,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(a_maxrownnz)*sizeof(**(&workj_L3)),(&workj_L3
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),791,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
792
793 /* This should be enough for almost all matrices. If not, memory is reallocated later. */
794 c_maxmem = 8*(ai[am]+bi[bm]);
795 /* Step 2: Populate pattern for C */
796 ierr = PetscMalloc1(c_maxmem,&cj)PetscMallocA(1,PETSC_FALSE,796,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(c_maxmem)*sizeof(**(&cj)),(&cj))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),796,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
797
798 ci_nnz = 0;
799 ci[0] = 0;
800 worki_L1[0] = 0;
801 worki_L2[0] = 0;
802 for (i=0; i<am; i++) {
803 const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
804 const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
805 rowsleft = anzi;
806 inputcol_L1 = acol;
807 L2_nnz = 0;
808 L2_nrows = 1; /* Number of rows to be merged on Level 3. output of L3 already exists -> initial value 1 */
809 worki_L2[1] = 0;
810 outputi_nnz = 0;
811
812 /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem -> allocate more memory */
813 while (ci_nnz+a_maxrownnz > c_maxmem) {
814 c_maxmem *= 2;
815 ndouble++;
816 ierr = PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj)((*PetscTrRealloc)((sizeof(PetscInt)*c_maxmem),816,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(void**)(&cj)))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),816,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
817 }
818
819 while (rowsleft) {
820 L1_rowsleft = PetscMin(64, rowsleft)(((64)<(rowsleft)) ? (64) : (rowsleft)); /* In the inner loop max 64 rows of B can be merged */
821 L1_nrows = 0;
822 L1_nnz = 0;
823 inputcol = inputcol_L1;
824 inputi = bi;
825 inputj = bj;
826
827 /* The following macro is used to specialize for small rows in A.
828 This helps with compiler unrolling, improving performance substantially.
829 Input: inputj inputi inputcol bn
830 Output: outputj outputi_nnz */
831 #define MatMatMultSymbolic_RowMergeMacro(ANNZ) \
832 window_min = bn; \
833 outputi_nnz = 0; \
834 for (k=0; k<ANNZ; ++k) { \
835 brow_ptr[k] = inputj + inputi[inputcol[k]]; \
836 brow_end[k] = inputj + inputi[inputcol[k]+1]; \
837 window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
838 window_min = PetscMin(window[k], window_min)(((window[k])<(window_min)) ? (window[k]) : (window_min)); \
839 } \
840 while (window_min < bn) { \
841 outputj[outputi_nnz++] = window_min; \
842 /* advance front and compute new minimum */ \
843 old_window_min = window_min; \
844 window_min = bn; \
845 for (k=0; k<ANNZ; ++k) { \
846 if (window[k] == old_window_min) { \
847 brow_ptr[k]++; \
848 window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
849 } \
850 window_min = PetscMin(window[k], window_min)(((window[k])<(window_min)) ? (window[k]) : (window_min)); \
851 } \
852 }
853
854 /************** L E V E L 1 ***************/
855 /* Merge up to 8 rows of B to L1 work array*/
856 while (L1_rowsleft) {
857 outputi_nnz = 0;
858 if (anzi > 8) outputj = workj_L1 + L1_nnz; /* Level 1 rowmerge*/
859 else outputj = cj + ci_nnz; /* Merge directly to C */
860
861 switch (L1_rowsleft) {
862 case 1: brow_ptr[0] = inputj + inputi[inputcol[0]];
863 brow_end[0] = inputj + inputi[inputcol[0]+1];
864 for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
865 inputcol += L1_rowsleft;
866 rowsleft -= L1_rowsleft;
867 L1_rowsleft = 0;
868 break;
869 case 2: MatMatMultSymbolic_RowMergeMacro(2);
870 inputcol += L1_rowsleft;
871 rowsleft -= L1_rowsleft;
872 L1_rowsleft = 0;
873 break;
874 case 3: MatMatMultSymbolic_RowMergeMacro(3);
875 inputcol += L1_rowsleft;
876 rowsleft -= L1_rowsleft;
877 L1_rowsleft = 0;
878 break;
879 case 4: MatMatMultSymbolic_RowMergeMacro(4);
880 inputcol += L1_rowsleft;
881 rowsleft -= L1_rowsleft;
882 L1_rowsleft = 0;
883 break;
884 case 5: MatMatMultSymbolic_RowMergeMacro(5);
885 inputcol += L1_rowsleft;
886 rowsleft -= L1_rowsleft;
887 L1_rowsleft = 0;
888 break;
889 case 6: MatMatMultSymbolic_RowMergeMacro(6);
890 inputcol += L1_rowsleft;
891 rowsleft -= L1_rowsleft;
892 L1_rowsleft = 0;
893 break;
894 case 7: MatMatMultSymbolic_RowMergeMacro(7);
895 inputcol += L1_rowsleft;
896 rowsleft -= L1_rowsleft;
897 L1_rowsleft = 0;
898 break;
899 default: MatMatMultSymbolic_RowMergeMacro(8);
900 inputcol += 8;
901 rowsleft -= 8;
902 L1_rowsleft -= 8;
903 break;
904 }
905 inputcol_L1 = inputcol;
906 L1_nnz += outputi_nnz;
907 worki_L1[++L1_nrows] = L1_nnz;
908 }
909
910 /********************** L E V E L 2 ************************/
911 /* Merge from L1 work array to either C or to L2 work array */
912 if (anzi > 8) {
913 inputi = worki_L1;
914 inputj = workj_L1;
915 inputcol = workcol;
916 outputi_nnz = 0;
917
918 if (anzi <= 64) outputj = cj + ci_nnz; /* Merge from L1 work array to C */
919 else outputj = workj_L2 + L2_nnz; /* Merge from L1 work array to L2 work array */
920
921 switch (L1_nrows) {
922 case 1: brow_ptr[0] = inputj + inputi[inputcol[0]];
923 brow_end[0] = inputj + inputi[inputcol[0]+1];
924 for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
925 break;
926 case 2: MatMatMultSymbolic_RowMergeMacro(2); break;
927 case 3: MatMatMultSymbolic_RowMergeMacro(3); break;
928 case 4: MatMatMultSymbolic_RowMergeMacro(4); break;
929 case 5: MatMatMultSymbolic_RowMergeMacro(5); break;
930 case 6: MatMatMultSymbolic_RowMergeMacro(6); break;
931 case 7: MatMatMultSymbolic_RowMergeMacro(7); break;
932 case 8: MatMatMultSymbolic_RowMergeMacro(8); break;
933 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatMult logic error: Not merging 1-8 rows from L1 work array!")return PetscError(((MPI_Comm)0x44000001),933,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,56,PETSC_ERROR_INITIAL,"MatMatMult logic error: Not merging 1-8 rows from L1 work array!"
)
;
934 }
935 L2_nnz += outputi_nnz;
936 worki_L2[++L2_nrows] = L2_nnz;
937
938 /************************ L E V E L 3 **********************/
939 /* Merge from L2 work array to either C or to L2 work array */
940 if (anzi > 64 && (L2_nrows == 8 || rowsleft == 0)) {
941 inputi = worki_L2;
942 inputj = workj_L2;
943 inputcol = workcol;
944 outputi_nnz = 0;
945 if (rowsleft) outputj = workj_L3;
946 else outputj = cj + ci_nnz;
947 switch (L2_nrows) {
948 case 1: brow_ptr[0] = inputj + inputi[inputcol[0]];
949 brow_end[0] = inputj + inputi[inputcol[0]+1];
950 for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
951 break;
952 case 2: MatMatMultSymbolic_RowMergeMacro(2); break;
953 case 3: MatMatMultSymbolic_RowMergeMacro(3); break;
954 case 4: MatMatMultSymbolic_RowMergeMacro(4); break;
955 case 5: MatMatMultSymbolic_RowMergeMacro(5); break;
956 case 6: MatMatMultSymbolic_RowMergeMacro(6); break;
957 case 7: MatMatMultSymbolic_RowMergeMacro(7); break;
958 case 8: MatMatMultSymbolic_RowMergeMacro(8); break;
959 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatMult logic error: Not merging 1-8 rows from L2 work array!")return PetscError(((MPI_Comm)0x44000001),959,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,56,PETSC_ERROR_INITIAL,"MatMatMult logic error: Not merging 1-8 rows from L2 work array!"
)
;
960 }
961 L2_nrows = 1;
962 L2_nnz = outputi_nnz;
963 worki_L2[1] = outputi_nnz;
964 /* Copy to workj_L2 */
965 if (rowsleft) {
966 for (k=0; k<outputi_nnz; ++k) workj_L2[k] = outputj[k];
967 }
968 }
969 }
970 } /* while (rowsleft) */
971#undef MatMatMultSymbolic_RowMergeMacro
972
973 /* terminate current row */
974 ci_nnz += outputi_nnz;
975 ci[i+1] = ci_nnz;
976 }
977
978 /* Step 3: Create the new symbolic matrix */
979 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),979,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
980 ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),980,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
981 ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),981,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
982
983 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
984 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
985 c = (Mat_SeqAIJ*)((*C)->data);
986 c->free_a = PETSC_TRUE;
987 c->free_ij = PETSC_TRUE;
988 c->nonew = 0;
989
990 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
991
992 /* set MatInfo */
993 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
994 if (afill < 1.0) afill = 1.0;
995 c->maxnz = ci[am];
996 c->nz = ci[am];
997 (*C)->info.mallocs = ndouble;
998 (*C)->info.fill_ratio_given = fill;
999 (*C)->info.fill_ratio_needed = afill;
1000
1001#if defined(PETSC_USE_INFO1)
1002 if (ci[am]) {
1003 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1003,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1004 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1004,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1005 } else {
1006 ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1006,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1007 }
1008#endif
1009
1010 /* Step 4: Free temporary work areas */
1011 ierr = PetscFree(workj_L1)((*PetscTrFree)((void*)(workj_L1),1011,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((workj_L1) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1011,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1012 ierr = PetscFree(workj_L2)((*PetscTrFree)((void*)(workj_L2),1012,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((workj_L2) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1012,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1013 ierr = PetscFree(workj_L3)((*PetscTrFree)((void*)(workj_L3),1013,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((workj_L3) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1013,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1014 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)
;
1015}
1016
1017/* concatenate unique entries and then sort */
1018PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,PetscReal fill,Mat *C)
1019{
1020 PetscErrorCode ierr;
1021 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
1022 const PetscInt *ai = a->i,*bi=b->i,*aj=a->j,*bj=b->j;
1023 PetscInt *ci,*cj;
1024 PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
1025 PetscReal afill;
1026 PetscInt i,j,ndouble = 0;
1027 PetscSegBuffer seg,segrow;
1028 char *seen;
1029
1030 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1030; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1031 ierr = PetscMalloc1(am+1,&ci)PetscMallocA(1,PETSC_FALSE,1031,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+1)*sizeof(**(&ci)),(&ci))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1031,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1032 ci[0] = 0;
1033
1034 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
1035 ierr = PetscSegBufferCreate(sizeof(PetscInt),(PetscInt)(fill*(ai[am]+bi[bm])),&seg);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1035,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1036 ierr = PetscSegBufferCreate(sizeof(PetscInt),100,&segrow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1036,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1037 ierr = PetscCalloc1(bn,&seen)PetscMallocA(1,PETSC_TRUE,1037,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(bn)*sizeof(**(&seen)),(&seen))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1037,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1038
1039 /* Determine ci and cj */
1040 for (i=0; i<am; i++) {
1041 const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
1042 const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
1043 PetscInt packlen = 0,*PETSC_RESTRICT__restrict crow;
1044 /* Pack segrow */
1045 for (j=0; j<anzi; j++) {
1046 PetscInt brow = acol[j],bjstart = bi[brow],bjend = bi[brow+1],k;
1047 for (k=bjstart; k<bjend; k++) {
1048 PetscInt bcol = bj[k];
1049 if (!seen[bcol]) { /* new entry */
1050 PetscInt *PETSC_RESTRICT__restrict slot;
1051 ierr = PetscSegBufferGetInts(segrow,1,&slot);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1051,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1052 *slot = bcol;
1053 seen[bcol] = 1;
1054 packlen++;
1055 }
1056 }
1057 }
1058 ierr = PetscSegBufferGetInts(seg,packlen,&crow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1058,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1059 ierr = PetscSegBufferExtractTo(segrow,crow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1059,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1060 ierr = PetscSortInt(packlen,crow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1060,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1061 ci[i+1] = ci[i] + packlen;
1062 for (j=0; j<packlen; j++) seen[crow[j]] = 0;
1063 }
1064 ierr = PetscSegBufferDestroy(&segrow);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1064,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1065 ierr = PetscFree(seen)((*PetscTrFree)((void*)(seen),1065,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((seen) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1065,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1066
1067 /* Column indices are in the segmented buffer */
1068 ierr = PetscSegBufferExtractAlloc(seg,&cj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1068,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1069 ierr = PetscSegBufferDestroy(&seg);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1069,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1070
1071 /* put together the new symbolic matrix */
1072 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1072,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1073 ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1073,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1074 ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1074,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1075
1076 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
1077 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
1078 c = (Mat_SeqAIJ*)((*C)->data);
1079 c->free_a = PETSC_TRUE;
1080 c->free_ij = PETSC_TRUE;
1081 c->nonew = 0;
1082
1083 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
1084
1085 /* set MatInfo */
1086 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
1087 if (afill < 1.0) afill = 1.0;
1088 c->maxnz = ci[am];
1089 c->nz = ci[am];
1090 (*C)->info.mallocs = ndouble;
1091 (*C)->info.fill_ratio_given = fill;
1092 (*C)->info.fill_ratio_needed = afill;
1093
1094#if defined(PETSC_USE_INFO1)
1095 if (ci[am]) {
1096 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill)PetscInfo_Private(__func__,(*C),"Reallocs %D; Fill ratio: given %g needed %g.\n"
,ndouble,(double)fill,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1096,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1097 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill)PetscInfo_Private(__func__,(*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n"
,(double)afill)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1097,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1098 } else {
1099 ierr = PetscInfo((*C),"Empty matrix product\n")PetscInfo_Private(__func__,(*C),"Empty matrix product\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1099,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1100 }
1101#endif
1102 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)
;
1103}
1104
1105/* This routine is not used. Should be removed! */
1106PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1107{
1108 PetscErrorCode ierr;
1109
1110 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1110; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1111 if (scall == MAT_INITIAL_MATRIX) {
1112 ierr = PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatTransposeMultSymbolic].active) ? (*PetscLogPLB)((MAT_MatTransposeMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1112,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1113 ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1113,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1114 ierr = PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatTransposeMultSymbolic].active) ? (*PetscLogPLE)((MAT_MatTransposeMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1114,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1115 }
1116 ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatTransposeMultNumeric].active) ? (*PetscLogPLB)((MAT_MatTransposeMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1116,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1117 ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1117,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1118 ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatTransposeMultNumeric].active) ? (*PetscLogPLE)((MAT_MatTransposeMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1118,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1119 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)
;
1120}
1121
1122PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(Mat A)
1123{
1124 PetscErrorCode ierr;
1125 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
1126 Mat_MatMatTransMult *abt=a->abt;
1127
1128 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1128; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1129 ierr = (abt->destroy)(A);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1129,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1130 ierr = MatTransposeColoringDestroy(&abt->matcoloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1130,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1131 ierr = MatDestroy(&abt->Bt_den);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1131,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1132 ierr = MatDestroy(&abt->ABt_den);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1132,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1133 ierr = PetscFree(abt)((*PetscTrFree)((void*)(abt),1133,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((abt) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1133,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1134 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)
;
1135}
1136
1137PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
1138{
1139 PetscErrorCode ierr;
1140 Mat Bt;
1141 PetscInt *bti,*btj;
1142 Mat_MatMatTransMult *abt;
1143 Mat_SeqAIJ *c;
1144
1145 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1145; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1146 /* create symbolic Bt */
1147 ierr = MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1147,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1148 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF((MPI_Comm)0x44000001),B->cmap->n,B->rmap->n,bti,btj,NULL((void*)0),&Bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1148,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1149 ierr = MatSetBlockSizes(Bt,PetscAbs(A->cmap->bs)(((A->cmap->bs) >= 0) ? (A->cmap->bs) : (-(A->
cmap->bs)))
,PetscAbs(B->cmap->bs)(((B->cmap->bs) >= 0) ? (B->cmap->bs) : (-(B->
cmap->bs)))
);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1149,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1150 ierr = MatSetType(Bt,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1150,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1151
1152 /* get symbolic C=A*Bt */
1153 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1153,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1154
1155 /* create a supporting struct for reuse intermidiate dense matrices with matcoloring */
1156 ierr = PetscNew(&abt)PetscMallocA(1,PETSC_TRUE,1156,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(1)*sizeof(**((&abt))),((&abt)))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1156,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1157 c = (Mat_SeqAIJ*)(*C)->data;
1158 c->abt = abt;
1159
1160 abt->usecoloring = PETSC_FALSE;
1161 abt->destroy = (*C)->ops->destroy;
1162 (*C)->ops->destroy = MatDestroy_SeqAIJ_MatMatMultTrans;
1163
1164 ierr = PetscOptionsGetBool(((PetscObject)A)->options,NULL((void*)0),"-matmattransmult_color",&abt->usecoloring,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1164,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1165 if (abt->usecoloring) {
1166 /* Create MatTransposeColoring from symbolic C=A*B^T */
1167 MatTransposeColoring matcoloring;
1168 MatColoring coloring;
1169 ISColoring iscoloring;
1170 Mat Bt_dense,C_dense;
1171 Mat_SeqAIJ *c=(Mat_SeqAIJ*)(*C)->data;
1172 /* inode causes memory problem, don't know why */
1173 if (c->inode.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MAT_USE_INODES is not supported. Use '-mat_no_inode'")return PetscError(((MPI_Comm)0x44000001),1173,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,56,PETSC_ERROR_INITIAL,"MAT_USE_INODES is not supported. Use '-mat_no_inode'"
)
;
1174
1175 ierr = MatColoringCreate(*C,&coloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1175,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1176 ierr = MatColoringSetDistance(coloring,2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1176,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1177 ierr = MatColoringSetType(coloring,MATCOLORINGSL"sl");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1177,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1178 ierr = MatColoringSetFromOptions(coloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1178,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1179 ierr = MatColoringApply(coloring,&iscoloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1179,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1180 ierr = MatColoringDestroy(&coloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1180,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1181 ierr = MatTransposeColoringCreate(*C,iscoloring,&matcoloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1181,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1182
1183 abt->matcoloring = matcoloring;
1184
1185 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1185,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1186
1187 /* Create Bt_dense and C_dense = A*Bt_dense */
1188 ierr = MatCreate(PETSC_COMM_SELF((MPI_Comm)0x44000001),&Bt_dense);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1188,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1189 ierr = MatSetSizes(Bt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1189,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1190 ierr = MatSetType(Bt_dense,MATSEQDENSE"seqdense");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1190,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1191 ierr = MatSeqDenseSetPreallocation(Bt_dense,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1191,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1192
1193 Bt_dense->assembled = PETSC_TRUE;
1194 abt->Bt_den = Bt_dense;
1195
1196 ierr = MatCreate(PETSC_COMM_SELF((MPI_Comm)0x44000001),&C_dense);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1196,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1197 ierr = MatSetSizes(C_dense,A->rmap->n,matcoloring->ncolors,A->rmap->n,matcoloring->ncolors);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1197,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1198 ierr = MatSetType(C_dense,MATSEQDENSE"seqdense");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1198,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1199 ierr = MatSeqDenseSetPreallocation(C_dense,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1199,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1200
1201 Bt_dense->assembled = PETSC_TRUE;
1202 abt->ABt_den = C_dense;
1203
1204#if defined(PETSC_USE_INFO1)
1205 {
1206 Mat_SeqAIJ *c = (Mat_SeqAIJ*)(*C)->data;
1207 ierr = PetscInfo7(*C,"Use coloring of C=A*B^T; B^T: %D %D, Bt_dense: %D,%D; Cnz %D / (cm*ncolors %D) = %g\n",B->cmap->n,B->rmap->n,Bt_dense->rmap->n,Bt_dense->cmap->n,c->nz,A->rmap->n*matcoloring->ncolors,(PetscReal)(c->nz)/(A->rmap->n*matcoloring->ncolors))PetscInfo_Private(__func__,*C,"Use coloring of C=A*B^T; B^T: %D %D, Bt_dense: %D,%D; Cnz %D / (cm*ncolors %D) = %g\n"
,B->cmap->n,B->rmap->n,Bt_dense->rmap->n,Bt_dense
->cmap->n,c->nz,A->rmap->n*matcoloring->ncolors
,(PetscReal)(c->nz)/(A->rmap->n*matcoloring->ncolors
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1207,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1208 }
1209#endif
1210 }
1211 /* clean up */
1212 ierr = MatDestroy(&Bt);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1212,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1213 ierr = MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1213,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1214 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)
;
1215}
1216
1217PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
1218{
1219 PetscErrorCode ierr;
1220 Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
1221 PetscInt *ai =a->i,*aj=a->j,*bi=b->i,*bj=b->j,anzi,bnzj,nexta,nextb,*acol,*bcol,brow;
1222 PetscInt cm =C->rmap->n,*ci=c->i,*cj=c->j,i,j,cnzi,*ccol;
1223 PetscLogDouble flops=0.0;
1224 MatScalar *aa =a->a,*aval,*ba=b->a,*bval,*ca,*cval;
1225 Mat_MatMatTransMult *abt = c->abt;
1226
1227 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1227; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1228 /* clear old values in C */
1229 if (!c->a) {
1230 ierr = PetscCalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_TRUE,1230,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[cm]+1)*sizeof(**(&ca)),(&ca))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1230,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1231 c->a = ca;
1232 c->free_a = PETSC_TRUE;
1233 } else {
1234 ca = c->a;
1235 ierr = PetscArrayzero(ca,ci[cm]+1)PetscMemzero(ca,(ci[cm]+1)*sizeof(*(ca)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1235,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1236 }
1237
1238 if (abt->usecoloring) {
1239 MatTransposeColoring matcoloring = abt->matcoloring;
1240 Mat Bt_dense,C_dense = abt->ABt_den;
1241
1242 /* Get Bt_dense by Apply MatTransposeColoring to B */
1243 Bt_dense = abt->Bt_den;
1244 ierr = MatTransColoringApplySpToDen(matcoloring,B,Bt_dense);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1244,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1245
1246 /* C_dense = A*Bt_dense */
1247 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,Bt_dense,C_dense);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1247,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1248
1249 /* Recover C from C_dense */
1250 ierr = MatTransColoringApplyDenToSp(matcoloring,C_dense,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1250,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1251 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)
;
1252 }
1253
1254 for (i=0; i<cm; i++) {
1255 anzi = ai[i+1] - ai[i];
1256 acol = aj + ai[i];
1257 aval = aa + ai[i];
1258 cnzi = ci[i+1] - ci[i];
1259 ccol = cj + ci[i];
1260 cval = ca + ci[i];
1261 for (j=0; j<cnzi; j++) {
1262 brow = ccol[j];
1263 bnzj = bi[brow+1] - bi[brow];
1264 bcol = bj + bi[brow];
1265 bval = ba + bi[brow];
1266
1267 /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */
1268 nexta = 0; nextb = 0;
1269 while (nexta<anzi && nextb<bnzj) {
1270 while (nexta < anzi && acol[nexta] < bcol[nextb]) nexta++;
1271 if (nexta == anzi) break;
1272 while (nextb < bnzj && acol[nexta] > bcol[nextb]) nextb++;
1273 if (nextb == bnzj) break;
1274 if (acol[nexta] == bcol[nextb]) {
1275 cval[j] += aval[nexta]*bval[nextb];
1276 nexta++; nextb++;
1277 flops += 2;
1278 }
1279 }
1280 }
1281 }
1282 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1282,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1283 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1283,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1284 ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1284,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1285 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)
;
1286}
1287
1288PetscErrorCode MatDestroy_SeqAIJ_MatTransMatMult(Mat A)
1289{
1290 PetscErrorCode ierr;
1291 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1292 Mat_MatTransMatMult *atb = a->atb;
1293
1294 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1294; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1295 if (atb) {
1296 ierr = MatDestroy(&atb->At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1296,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1297 ierr = (*atb->destroy)(A);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1297,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1298 }
1299 ierr = PetscFree(atb)((*PetscTrFree)((void*)(atb),1299,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((atb) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1299,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1300 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)
;
1301}
1302
1303PetscErrorCode MatTransposeMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1304{
1305 PetscErrorCode ierr;
1306 const char *algTypes[2] = {"matmatmult","outerproduct"};
1307 PetscInt alg=0; /* set default algorithm */
1308 Mat At;
1309 Mat_MatTransMatMult *atb;
1310 Mat_SeqAIJ *c;
1311
1312 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1312; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1313 if (scall == MAT_INITIAL_MATRIX) {
1314 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatTransposeMatMult","Mat")0; do { PetscOptionItems PetscOptionsObjectBase; PetscOptionItems
*PetscOptionsObject = &PetscOptionsObjectBase; PetscMemzero
(PetscOptionsObject,sizeof(PetscOptionItems)); for (PetscOptionsObject
->count=(PetscOptionsPublish?-1:1); PetscOptionsObject->
count<2; PetscOptionsObject->count++) { PetscErrorCode _5_ierr
= PetscOptionsBegin_Private(PetscOptionsObject,PetscObjectComm
((PetscObject)A),((PetscObject)A)->prefix,"MatTransposeMatMult"
,"Mat");do {if (__builtin_expect(!!(_5_ierr),0)) return PetscError
(((MPI_Comm)0x44000001),1314,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1314,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1315 ierr = PetscOptionsEList("-mattransposematmult_via","Algorithmic approach","MatTransposeMatMult",algTypes,2,algTypes[0],&alg,NULL)PetscOptionsEList_Private(PetscOptionsObject,"-mattransposematmult_via"
,"Algorithmic approach","MatTransposeMatMult",algTypes,2,algTypes
[0],&alg,((void*)0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1315,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1316 ierr = PetscOptionsEnd()_5_ierr = PetscOptionsEnd_Private(PetscOptionsObject);do {if (
__builtin_expect(!!(_5_ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1316,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);}} while (0)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1316,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1317
1318 switch (alg) {
1319 case 1:
1320 ierr = MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1320,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1321 break;
1322 default:
1323 ierr = PetscNew(&atb)PetscMallocA(1,PETSC_TRUE,1323,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(1)*sizeof(**((&atb))),((&atb)))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1323,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1324 ierr = MatTranspose_SeqAIJ(A,MAT_INITIAL_MATRIX,&At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1324,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1325 ierr = MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_INITIAL_MATRIX,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1325,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1326
1327 c = (Mat_SeqAIJ*)(*C)->data;
1328 c->atb = atb;
1329 atb->At = At;
1330 atb->destroy = (*C)->ops->destroy;
1331 (*C)->ops->destroy = MatDestroy_SeqAIJ_MatTransMatMult;
1332
1333 break;
1334 }
1335 }
1336 if (alg) {
1337 ierr = (*(*C)->ops->mattransposemultnumeric)(A,B,*C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1337,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1338 } else if (!alg && scall == MAT_REUSE_MATRIX) {
1339 c = (Mat_SeqAIJ*)(*C)->data;
1340 atb = c->atb;
1341 At = atb->At;
1342 ierr = MatTranspose_SeqAIJ(A,MAT_REUSE_MATRIX,&At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1342,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1343 ierr = MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_REUSE_MATRIX,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1343,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1344 }
1345 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)
;
1346}
1347
1348PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
1349{
1350 PetscErrorCode ierr;
1351 Mat At;
1352 PetscInt *ati,*atj;
1353
1354 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1354; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1355 /* create symbolic At */
1356 ierr = MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1356,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1357 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF((MPI_Comm)0x44000001),A->cmap->n,A->rmap->n,ati,atj,NULL((void*)0),&At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1357,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1358 ierr = MatSetBlockSizes(At,PetscAbs(A->cmap->bs)(((A->cmap->bs) >= 0) ? (A->cmap->bs) : (-(A->
cmap->bs)))
,PetscAbs(B->cmap->bs)(((B->cmap->bs) >= 0) ? (B->cmap->bs) : (-(B->
cmap->bs)))
);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1358,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1359 ierr = MatSetType(At,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1359,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1360
1361 /* get symbolic C=At*B */
1362 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1362,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1363
1364 /* clean up */
1365 ierr = MatDestroy(&At);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1365,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1366 ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1366,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1367
1368 (*C)->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ;
1369 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)
;
1370}
1371
1372PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
1373{
1374 PetscErrorCode ierr;
1375 Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
1376 PetscInt am =A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb;
1377 PetscInt cm =C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k;
1378 PetscLogDouble flops=0.0;
1379 MatScalar *aa =a->a,*ba,*ca,*caj;
1380
1381 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1381; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1382 if (!c->a) {
1383 ierr = PetscCalloc1(ci[cm]+1,&ca)PetscMallocA(1,PETSC_TRUE,1383,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(ci[cm]+1)*sizeof(**(&ca)),(&ca))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1383,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1384
1385 c->a = ca;
1386 c->free_a = PETSC_TRUE;
1387 } else {
1388 ca = c->a;
1389 ierr = PetscArrayzero(ca,ci[cm])PetscMemzero(ca,(ci[cm])*sizeof(*(ca)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1389,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1390 }
1391
1392 /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */
1393 for (i=0; i<am; i++) {
1394 bj = b->j + bi[i];
1395 ba = b->a + bi[i];
1396 bnzi = bi[i+1] - bi[i];
1397 anzi = ai[i+1] - ai[i];
1398 for (j=0; j<anzi; j++) {
1399 nextb = 0;
1400 crow = *aj++;
1401 cjj = cj + ci[crow];
1402 caj = ca + ci[crow];
1403 /* perform sparse axpy operation. Note cjj includes bj. */
1404 for (k=0; nextb<bnzi; k++) {
1405 if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */
1406 caj[k] += (*aa)*(*(ba+nextb));
1407 nextb++;
1408 }
1409 }
1410 flops += 2*bnzi;
1411 aa++;
1412 }
1413 }
1414
1415 /* Assemble the final matrix and clean up */
1416 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1416,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1417 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1417,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1418 ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1418,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1419 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)
;
1420}
1421
1422PETSC_INTERNextern __attribute__((visibility ("hidden"))) PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1423{
1424 PetscErrorCode ierr;
1425
1426 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1426; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1427 if (scall == MAT_INITIAL_MATRIX) {
1428 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultSymbolic].active) ? (*PetscLogPLB)((MAT_MatMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1428,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1429 ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1429,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1430 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultSymbolic].active) ? (*PetscLogPLE)((MAT_MatMultSymbolic
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1430,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1431 }
1432 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0)(((PetscLogPLB && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultNumeric].active) ? (*PetscLogPLB)((MAT_MatMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1432,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1433 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1433,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1434 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0)(((PetscLogPLE && petsc_stageLog->stageInfo[petsc_stageLog
->curStage].perfInfo.active && petsc_stageLog->
stageInfo[petsc_stageLog->curStage].eventLog->eventInfo
[MAT_MatMultNumeric].active) ? (*PetscLogPLE)((MAT_MatMultNumeric
),0,(PetscObject)(A),(PetscObject)(B),(PetscObject)(0),(PetscObject
)(0)) : 0 ))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1434,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1435 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)
;
1436}
1437
1438PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1439{
1440 PetscErrorCode ierr;
1441
1442 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1442; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1443 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1443,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1444
1445 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1446 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)
;
1447}
1448
1449PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
1450{
1451 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
1452 Mat_SeqDense *bd = (Mat_SeqDense*)B->data;
1453 PetscErrorCode ierr;
1454 PetscScalar *c,*b,r1,r2,r3,r4,*c1,*c2,*c3,*c4,aatmp;
1455 const PetscScalar *aa,*b1,*b2,*b3,*b4;
1456 const PetscInt *aj;
1457 PetscInt cm=C->rmap->n,cn=B->cmap->n,bm=bd->lda,am=A->rmap->n;
1458 PetscInt am4=4*am,bm4=4*bm,col,i,j,n,ajtmp;
1459
1460 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1460; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1461 if (!cm || !cn) PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0)
;
1462 if (B->rmap->n != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in A %D not equal rows in B %D\n",A->cmap->n,B->rmap->n)return PetscError(((MPI_Comm)0x44000001),1462,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,60,PETSC_ERROR_INITIAL,"Number columns in A %D not equal rows in B %D\n"
,A->cmap->n,B->rmap->n)
;
1463 if (A->rmap->n != C->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows in C %D not equal rows in A %D\n",C->rmap->n,A->rmap->n)return PetscError(((MPI_Comm)0x44000001),1463,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,60,PETSC_ERROR_INITIAL,"Number rows in C %D not equal rows in A %D\n"
,C->rmap->n,A->rmap->n)
;
1464 if (B->cmap->n != C->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in B %D not equal columns in C %D\n",B->cmap->n,C->cmap->n)return PetscError(((MPI_Comm)0x44000001),1464,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,60,PETSC_ERROR_INITIAL,"Number columns in B %D not equal columns in C %D\n"
,B->cmap->n,C->cmap->n)
;
1465 b = bd->v;
1466 ierr = MatDenseGetArray(C,&c);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1466,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1467 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
1468 c1 = c; c2 = c1 + am; c3 = c2 + am; c4 = c3 + am;
1469 for (col=0; col<cn-4; col += 4) { /* over columns of C */
1470 for (i=0; i<am; i++) { /* over rows of C in those columns */
1471 r1 = r2 = r3 = r4 = 0.0;
1472 n = a->i[i+1] - a->i[i];
1473 aj = a->j + a->i[i];
1474 aa = a->a + a->i[i];
1475 for (j=0; j<n; j++) {
1476 aatmp = aa[j]; ajtmp = aj[j];
1477 r1 += aatmp*b1[ajtmp];
1478 r2 += aatmp*b2[ajtmp];
1479 r3 += aatmp*b3[ajtmp];
1480 r4 += aatmp*b4[ajtmp];
1481 }
1482 c1[i] = r1;
1483 c2[i] = r2;
1484 c3[i] = r3;
1485 c4[i] = r4;
1486 }
1487 b1 += bm4; b2 += bm4; b3 += bm4; b4 += bm4;
1488 c1 += am4; c2 += am4; c3 += am4; c4 += am4;
1489 }
1490 for (; col<cn; col++) { /* over extra columns of C */
1491 for (i=0; i<am; i++) { /* over rows of C in those columns */
1492 r1 = 0.0;
1493 n = a->i[i+1] - a->i[i];
1494 aj = a->j + a->i[i];
1495 aa = a->a + a->i[i];
1496 for (j=0; j<n; j++) {
1497 r1 += aa[j]*b1[aj[j]];
1498 }
1499 c1[i] = r1;
1500 }
1501 b1 += bm;
1502 c1 += am;
1503 }
1504 ierr = PetscLogFlops(cn*(2.0*a->nz));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1504,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1505 ierr = MatDenseRestoreArray(C,&c);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1505,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1506 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1506,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1507 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1507,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1508 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)
;
1509}
1510
1511/*
1512 Note very similar to MatMult_SeqAIJ(), should generate both codes from same base
1513*/
1514PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
1515{
1516 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
1517 Mat_SeqDense *bd = (Mat_SeqDense*)B->data;
1518 PetscErrorCode ierr;
1519 PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
1520 MatScalar *aa;
1521 PetscInt cm = C->rmap->n, cn=B->cmap->n, bm=bd->lda, col, i,j,n,*aj, am = A->rmap->n,*ii,arm;
1522 PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam,*ridx;
1523
1524 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1524; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1525 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)
;
1526 b = bd->v;
1527 ierr = MatDenseGetArray(C,&c);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1527,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1528 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
1529
1530 if (a->compressedrow.use) { /* use compressed row format */
1531 for (col=0; col<cn-4; col += 4) { /* over columns of C */
1532 colam = col*am;
1533 arm = a->compressedrow.nrows;
1534 ii = a->compressedrow.i;
1535 ridx = a->compressedrow.rindex;
1536 for (i=0; i<arm; i++) { /* over rows of C in those columns */
1537 r1 = r2 = r3 = r4 = 0.0;
1538 n = ii[i+1] - ii[i];
1539 aj = a->j + ii[i];
1540 aa = a->a + ii[i];
1541 for (j=0; j<n; j++) {
1542 r1 += (*aa)*b1[*aj];
1543 r2 += (*aa)*b2[*aj];
1544 r3 += (*aa)*b3[*aj];
1545 r4 += (*aa++)*b4[*aj++];
1546 }
1547 c[colam + ridx[i]] += r1;
1548 c[colam + am + ridx[i]] += r2;
1549 c[colam + am2 + ridx[i]] += r3;
1550 c[colam + am3 + ridx[i]] += r4;
1551 }
1552 b1 += bm4;
1553 b2 += bm4;
1554 b3 += bm4;
1555 b4 += bm4;
1556 }
1557 for (; col<cn; col++) { /* over extra columns of C */
1558 colam = col*am;
1559 arm = a->compressedrow.nrows;
1560 ii = a->compressedrow.i;
1561 ridx = a->compressedrow.rindex;
1562 for (i=0; i<arm; i++) { /* over rows of C in those columns */
1563 r1 = 0.0;
1564 n = ii[i+1] - ii[i];
1565 aj = a->j + ii[i];
1566 aa = a->a + ii[i];
1567
1568 for (j=0; j<n; j++) {
1569 r1 += (*aa++)*b1[*aj++];
1570 }
1571 c[colam + ridx[i]] += r1;
1572 }
1573 b1 += bm;
1574 }
1575 } else {
1576 for (col=0; col<cn-4; col += 4) { /* over columns of C */
1577 colam = col*am;
1578 for (i=0; i<am; i++) { /* over rows of C in those columns */
1579 r1 = r2 = r3 = r4 = 0.0;
1580 n = a->i[i+1] - a->i[i];
1581 aj = a->j + a->i[i];
1582 aa = a->a + a->i[i];
1583 for (j=0; j<n; j++) {
1584 r1 += (*aa)*b1[*aj];
1585 r2 += (*aa)*b2[*aj];
1586 r3 += (*aa)*b3[*aj];
1587 r4 += (*aa++)*b4[*aj++];
1588 }
1589 c[colam + i] += r1;
1590 c[colam + am + i] += r2;
1591 c[colam + am2 + i] += r3;
1592 c[colam + am3 + i] += r4;
1593 }
1594 b1 += bm4;
1595 b2 += bm4;
1596 b3 += bm4;
1597 b4 += bm4;
1598 }
1599 for (; col<cn; col++) { /* over extra columns of C */
1600 colam = col*am;
1601 for (i=0; i<am; i++) { /* over rows of C in those columns */
1602 r1 = 0.0;
1603 n = a->i[i+1] - a->i[i];
1604 aj = a->j + a->i[i];
1605 aa = a->a + a->i[i];
1606
1607 for (j=0; j<n; j++) {
1608 r1 += (*aa++)*b1[*aj++];
1609 }
1610 c[colam + i] += r1;
1611 }
1612 b1 += bm;
1613 }
1614 }
1615 ierr = PetscLogFlops(cn*2.0*a->nz);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1615,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1616 ierr = MatDenseRestoreArray(C,&c);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1616,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1617 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)
;
1618}
1619
1620PetscErrorCode MatTransColoringApplySpToDen_SeqAIJ(MatTransposeColoring coloring,Mat B,Mat Btdense)
1621{
1622 PetscErrorCode ierr;
1623 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
1624 Mat_SeqDense *btdense = (Mat_SeqDense*)Btdense->data;
1625 PetscInt *bi = b->i,*bj=b->j;
1626 PetscInt m = Btdense->rmap->n,n=Btdense->cmap->n,j,k,l,col,anz,*btcol,brow,ncolumns;
1627 MatScalar *btval,*btval_den,*ba=b->a;
1628 PetscInt *columns=coloring->columns,*colorforcol=coloring->colorforcol,ncolors=coloring->ncolors;
1629
1630 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1630; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1631 btval_den=btdense->v;
1632 ierr = PetscArrayzero(btval_den,m*n)PetscMemzero(btval_den,(m*n)*sizeof(*(btval_den)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1632,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1633 for (k=0; k<ncolors; k++) {
1634 ncolumns = coloring->ncolumns[k];
1635 for (l=0; l<ncolumns; l++) { /* insert a row of B to a column of Btdense */
1636 col = *(columns + colorforcol[k] + l);
1637 btcol = bj + bi[col];
1638 btval = ba + bi[col];
1639 anz = bi[col+1] - bi[col];
1640 for (j=0; j<anz; j++) {
1641 brow = btcol[j];
1642 btval_den[brow] = btval[j];
1643 }
1644 }
1645 btval_den += m;
1646 }
1647 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)
;
1648}
1649
1650PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring,Mat Cden,Mat Csp)
1651{
1652 PetscErrorCode ierr;
1653 Mat_SeqAIJ *csp = (Mat_SeqAIJ*)Csp->data;
1654 const PetscScalar *ca_den,*ca_den_ptr;
1655 PetscScalar *ca=csp->a;
1656 PetscInt k,l,m=Cden->rmap->n,ncolors=matcoloring->ncolors;
1657 PetscInt brows=matcoloring->brows,*den2sp=matcoloring->den2sp;
1658 PetscInt nrows,*row,*idx;
1659 PetscInt *rows=matcoloring->rows,*colorforrow=matcoloring->colorforrow;
1660
1661 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1661; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1662 ierr = MatDenseGetArrayRead(Cden,&ca_den);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1662,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1663
1664 if (brows > 0) {
1665 PetscInt *lstart,row_end,row_start;
1666 lstart = matcoloring->lstart;
1667 ierr = PetscArrayzero(lstart,ncolors)PetscMemzero(lstart,(ncolors)*sizeof(*(lstart)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1667,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1668
1669 row_end = brows;
1670 if (row_end > m) row_end = m;
1671 for (row_start=0; row_start<m; row_start+=brows) { /* loop over row blocks of Csp */
1672 ca_den_ptr = ca_den;
1673 for (k=0; k<ncolors; k++) { /* loop over colors (columns of Cden) */
1674 nrows = matcoloring->nrows[k];
1675 row = rows + colorforrow[k];
1676 idx = den2sp + colorforrow[k];
1677 for (l=lstart[k]; l<nrows; l++) {
1678 if (row[l] >= row_end) {
1679 lstart[k] = l;
1680 break;
1681 } else {
1682 ca[idx[l]] = ca_den_ptr[row[l]];
1683 }
1684 }
1685 ca_den_ptr += m;
1686 }
1687 row_end += brows;
1688 if (row_end > m) row_end = m;
1689 }
1690 } else { /* non-blocked impl: loop over columns of Csp - slow if Csp is large */
1691 ca_den_ptr = ca_den;
1692 for (k=0; k<ncolors; k++) {
1693 nrows = matcoloring->nrows[k];
1694 row = rows + colorforrow[k];
1695 idx = den2sp + colorforrow[k];
1696 for (l=0; l<nrows; l++) {
1697 ca[idx[l]] = ca_den_ptr[row[l]];
1698 }
1699 ca_den_ptr += m;
1700 }
1701 }
1702
1703 ierr = MatDenseRestoreArrayRead(Cden,&ca_den);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1703,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1704#if defined(PETSC_USE_INFO1)
1705 if (matcoloring->brows > 0) {
1706 ierr = PetscInfo1(Csp,"Loop over %D row blocks for den2sp\n",brows)PetscInfo_Private(__func__,Csp,"Loop over %D row blocks for den2sp\n"
,brows)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1706,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1707 } else {
1708 ierr = PetscInfo(Csp,"Loop over colors/columns of Cden, inefficient for large sparse matrix product \n")PetscInfo_Private(__func__,Csp,"Loop over colors/columns of Cden, inefficient for large sparse matrix product \n"
)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1708,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1709 }
1710#endif
1711 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)
;
1712}
1713
1714PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatTransposeColoring c)
1715{
1716 PetscErrorCode ierr;
1717 PetscInt i,n,nrows,Nbs,j,k,m,ncols,col,cm;
1718 const PetscInt *is,*ci,*cj,*row_idx;
1719 PetscInt nis = iscoloring->n,*rowhit,bs = 1;
1720 IS *isa;
1721 Mat_SeqAIJ *csp = (Mat_SeqAIJ*)mat->data;
1722 PetscInt *colorforrow,*rows,*rows_i,*idxhit,*spidx,*den2sp,*den2sp_i;
1723 PetscInt *colorforcol,*columns,*columns_i,brows;
1724 PetscBool flg;
1725
1726 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1726; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1727 ierr = ISColoringGetIS(iscoloring,PETSC_USE_POINTER,PETSC_IGNORE((void*)0),&isa);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1727,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1728
1729 /* bs >1 is not being tested yet! */
1730 Nbs = mat->cmap->N/bs;
1731 c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */
1732 c->N = Nbs;
1733 c->m = c->M;
1734 c->rstart = 0;
1735 c->brows = 100;
1736
1737 c->ncolors = nis;
1738 ierr = PetscMalloc3(nis,&c->ncolumns,nis,&c->nrows,nis+1,&colorforrow)PetscMallocA(3,PETSC_FALSE,1738,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(nis)*sizeof(**(&c->ncolumns)),(&c->ncolumns
),(size_t)(nis)*sizeof(**(&c->nrows)),(&c->nrows
),(size_t)(nis+1)*sizeof(**(&colorforrow)),(&colorforrow
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1738,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1739 ierr = PetscMalloc1(csp->nz+1,&rows)PetscMallocA(1,PETSC_FALSE,1739,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(csp->nz+1)*sizeof(**(&rows)),(&rows))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1739,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1740 ierr = PetscMalloc1(csp->nz+1,&den2sp)PetscMallocA(1,PETSC_FALSE,1740,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(csp->nz+1)*sizeof(**(&den2sp)),(&den2sp)
)
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1740,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1741
1742 brows = c->brows;
1743 ierr = PetscOptionsGetInt(NULL((void*)0),NULL((void*)0),"-matden2sp_brows",&brows,&flg);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1743,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1744 if (flg) c->brows = brows;
1745 if (brows > 0) {
1746 ierr = PetscMalloc1(nis+1,&c->lstart)PetscMallocA(1,PETSC_FALSE,1746,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(nis+1)*sizeof(**(&c->lstart)),(&c->lstart
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1746,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1747 }
1748
1749 colorforrow[0] = 0;
1750 rows_i = rows;
1751 den2sp_i = den2sp;
1752
1753 ierr = PetscMalloc1(nis+1,&colorforcol)PetscMallocA(1,PETSC_FALSE,1753,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(nis+1)*sizeof(**(&colorforcol)),(&colorforcol
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1753,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1754 ierr = PetscMalloc1(Nbs+1,&columns)PetscMallocA(1,PETSC_FALSE,1754,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(Nbs+1)*sizeof(**(&columns)),(&columns))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1754,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1755
1756 colorforcol[0] = 0;
1757 columns_i = columns;
1758
1759 /* get column-wise storage of mat */
1760 ierr = MatGetColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1760,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1761
1762 cm = c->m;
1763 ierr = PetscMalloc1(cm+1,&rowhit)PetscMallocA(1,PETSC_FALSE,1763,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(cm+1)*sizeof(**(&rowhit)),(&rowhit))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1763,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1764 ierr = PetscMalloc1(cm+1,&idxhit)PetscMallocA(1,PETSC_FALSE,1764,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(cm+1)*sizeof(**(&idxhit)),(&idxhit))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1764,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1765 for (i=0; i<nis; i++) { /* loop over color */
1766 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1766,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1767 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1767,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1768
1769 c->ncolumns[i] = n;
1770 if (n) {
1771 ierr = PetscArraycpy(columns_i,is,n)((sizeof(*(columns_i)) != sizeof(*(is))) || PetscMemcpy(columns_i
,is,(n)*sizeof(*(columns_i))))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1771,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1772 }
1773 colorforcol[i+1] = colorforcol[i] + n;
1774 columns_i += n;
1775
1776 /* fast, crude version requires O(N*N) work */
1777 ierr = PetscArrayzero(rowhit,cm)PetscMemzero(rowhit,(cm)*sizeof(*(rowhit)));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1777,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1778
1779 for (j=0; j<n; j++) { /* loop over columns*/
1780 col = is[j];
1781 row_idx = cj + ci[col];
1782 m = ci[col+1] - ci[col];
1783 for (k=0; k<m; k++) { /* loop over columns marking them in rowhit */
1784 idxhit[*row_idx] = spidx[ci[col] + k];
1785 rowhit[*row_idx++] = col + 1;
1786 }
1787 }
1788 /* count the number of hits */
1789 nrows = 0;
1790 for (j=0; j<cm; j++) {
1791 if (rowhit[j]) nrows++;
1792 }
1793 c->nrows[i] = nrows;
1794 colorforrow[i+1] = colorforrow[i] + nrows;
1795
1796 nrows = 0;
1797 for (j=0; j<cm; j++) { /* loop over rows */
1798 if (rowhit[j]) {
1799 rows_i[nrows] = j;
1800 den2sp_i[nrows] = idxhit[j];
1801 nrows++;
1802 }
1803 }
1804 den2sp_i += nrows;
1805
1806 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1806,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1807 rows_i += nrows;
1808 }
1809 ierr = MatRestoreColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1809,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1810 ierr = PetscFree(rowhit)((*PetscTrFree)((void*)(rowhit),1810,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((rowhit) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1810,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1811 ierr = ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&isa);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1811,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1812 if (csp->nz != colorforrow[nis]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"csp->nz %d != colorforrow[nis] %d",csp->nz,colorforrow[nis])return PetscError(((MPI_Comm)0x44000001),1812,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,77,PETSC_ERROR_INITIAL,"csp->nz %d != colorforrow[nis] %d"
,csp->nz,colorforrow[nis])
;
1813
1814 c->colorforrow = colorforrow;
1815 c->rows = rows;
1816 c->den2sp = den2sp;
1817 c->colorforcol = colorforcol;
1818 c->columns = columns;
1819
1820 ierr = PetscFree(idxhit)((*PetscTrFree)((void*)(idxhit),1820,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((idxhit) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1820,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1821 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)
;
1822}
1823
1824/* The combine method has done the symbolic and numeric in the first phase, and so we just return here */
1825PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Combined(Mat A,Mat B,Mat C)
1826{
1827 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1827; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1828 /* Empty function */
1829 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)
;
1830}
1831
1832/* This algorithm combines the symbolic and numeric phase of matrix-matrix multiplication. */
1833PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ_Combined(Mat A,Mat B,PetscReal fill,Mat *C)
1834{
1835 PetscErrorCode ierr;
1836 PetscLogDouble flops=0.0;
1837 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
1838 const PetscInt *ai=a->i,*bi=b->i;
1839 PetscInt *ci,*cj,*cj_i;
1840 PetscScalar *ca,*ca_i;
1841 PetscInt b_maxmemrow,c_maxmem,a_col;
1842 PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
1843 PetscInt i,k,ndouble=0;
1844 PetscReal afill;
1845 PetscScalar *c_row_val_dense;
1846 PetscBool *c_row_idx_flags;
1847 PetscInt *aj_i=a->j;
1848 PetscScalar *aa_i=a->a;
1849
1850 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
; petscstack->line[petscstack->currentsize] = 1850; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0)
;
1851
1852 /* Step 1: Determine upper bounds on memory for C and allocate memory */
1853 /* This should be enough for almost all matrices. If still more memory is needed, it is reallocated later. */
1854 c_maxmem = 8*(ai[am]+bi[bm]);
1855 b_maxmemrow = PetscMin(bi[bm],bn)(((bi[bm])<(bn)) ? (bi[bm]) : (bn));
1856 ierr = PetscMalloc1(am+1,&ci)PetscMallocA(1,PETSC_FALSE,1856,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(am+1)*sizeof(**(&ci)),(&ci))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1856,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1857 ierr = PetscCalloc1(bn,&c_row_val_dense)PetscMallocA(1,PETSC_TRUE,1857,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(bn)*sizeof(**(&c_row_val_dense)),(&c_row_val_dense
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1857,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1858 ierr = PetscCalloc1(bn,&c_row_idx_flags)PetscMallocA(1,PETSC_TRUE,1858,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(bn)*sizeof(**(&c_row_idx_flags)),(&c_row_idx_flags
))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1858,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1859 ierr = PetscMalloc1(c_maxmem,&cj)PetscMallocA(1,PETSC_FALSE,1859,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(c_maxmem)*sizeof(**(&cj)),(&cj))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1859,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1860 ierr = PetscMalloc1(c_maxmem,&ca)PetscMallocA(1,PETSC_FALSE,1860,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(size_t)(c_maxmem)*sizeof(**(&ca)),(&ca))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1860,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1861 ca_i = ca;
1862 cj_i = cj;
1863 ci[0] = 0;
1864 for (i=0; i<am; i++) {
1865 /* Step 2: Initialize the dense row vector for C */
1866 const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
1867 PetscInt cnzi = 0;
1868 PetscInt *bj_i;
1869 PetscScalar *ba_i;
1870 /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem -> allocate more memory
1871 Usually, there is enough memory in the first place, so this is not executed. */
1872 while (ci[i] + b_maxmemrow > c_maxmem) {
1873 c_maxmem *= 2;
1874 ndouble++;
1875 ierr = PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj)((*PetscTrRealloc)((sizeof(PetscInt)*c_maxmem),1875,__func__,
"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(void**)(&cj)))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1875,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1876 ierr = PetscRealloc(sizeof(PetscScalar)*c_maxmem,&ca)((*PetscTrRealloc)((sizeof(PetscScalar)*c_maxmem),1876,__func__
,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,(void**)(&ca)))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1876,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1877 }
1878
1879 /* Step 3: Do the numerical calculations */
1880 for (a_col=0; a_col<anzi; a_col++) { /* iterate over all non zero values in a row of A */
1881 PetscInt a_col_index = aj_i[a_col];
1882 const PetscInt bnzi = bi[a_col_index+1] - bi[a_col_index];
1883 flops += 2*bnzi;
1884 bj_i = b->j + bi[a_col_index]; /* points to the current row in bj */
1885 ba_i = b->a + bi[a_col_index]; /* points to the current row in ba */
1886 for (k=0; k<bnzi; ++k) { /* iterate over all non zeros of this row in B */
1887 if (c_row_idx_flags[bj_i[k]] == PETSC_FALSE) {
1888 cj_i[cnzi++] = bj_i[k];
1889 c_row_idx_flags[bj_i[k]] = PETSC_TRUE;
1890 }
1891 c_row_val_dense[bj_i[k]] += aa_i[a_col] * ba_i[k];
1892 }
1893 }
1894
1895 /* Sort array */
1896 ierr = PetscSortInt(cnzi,cj_i);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1896,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1897 /* Step 4 */
1898 for (k=0; k<cnzi; k++) {
1899 ca_i[k] = c_row_val_dense[cj_i[k]];
1900 c_row_val_dense[cj_i[k]] = 0.;
1901 c_row_idx_flags[cj_i[k]] = PETSC_FALSE;
1902 }
1903 /* terminate current row */
1904 aa_i += anzi;
1905 aj_i += anzi;
1906 ca_i += cnzi;
1907 cj_i += cnzi;
1908 ci[i+1] = ci[i] + cnzi;
1909 flops += cnzi;
1910 }
1911
1912 /* Step 5 */
1913 /* Create the new matrix */
1914 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL((void*)0),C);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1914,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1915 ierr = MatSetBlockSizesFromMats(*C,A,B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1915,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1916 ierr = MatSetType(*C,((PetscObject)A)->type_name);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1916,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1917
1918 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
1919 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
1920 c = (Mat_SeqAIJ*)((*C)->data);
1921 c->a = ca;
1922 c->free_a = PETSC_TRUE;
1923 c->free_ij = PETSC_TRUE;
1924 c->nonew = 0;
1925
1926 /* set MatInfo */
1927 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
1928 if (afill < 1.0) afill = 1.0;
1929 c->maxnz = ci[am];
1930 c->nz = ci[am];
1931 (*C)->info.mallocs = ndouble;
1932 (*C)->info.fill_ratio_given = fill;
1933 (*C)->info.fill_ratio_needed = afill;
1934 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Combined;
1935
1936 ierr = PetscFree(c_row_val_dense)((*PetscTrFree)((void*)(c_row_val_dense),1936,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((c_row_val_dense) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1936,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1937 ierr = PetscFree(c_row_idx_flags)((*PetscTrFree)((void*)(c_row_idx_flags),1937,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
) || ((c_row_idx_flags) = 0,0))
;CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1937,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1938
1939 ierr = MatAssemblyBegin(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1939,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1940 ierr = MatAssemblyEnd(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1940,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1941 ierr = PetscLogFlops(flops);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),1941,__func__,"/sandbox/petsc/petsc.next-tmp/src/mat/impls/aij/seq/matmatmult.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0)
;
1942 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)
;
1943}