1 #define PETSCMAT_DLL 2 3 /* 4 Defines matrix-matrix product routines for pairs of SeqAIJ matrices 5 C = A * B 6 */ 7 8 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 9 #include "src/mat/utils/freespace.h" 10 #include "petscbt.h" 11 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqAIJ" 15 PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 16 { 17 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 if (scall == MAT_INITIAL_MATRIX){ 21 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 22 } 23 ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 24 PetscFunctionReturn(0); 25 } 26 27 28 #undef __FUNCT__ 29 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqAIJ" 30 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 31 { 32 PetscErrorCode ierr; 33 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 34 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 35 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci,*cj; 36 PetscInt am=A->M,bn=B->N,bm=B->M; 37 PetscInt i,j,anzi,brow,bnzj,cnzi,nlnk,*lnk,nspacedouble=0; 38 MatScalar *ca; 39 PetscBT lnkbt; 40 41 PetscFunctionBegin; 42 /* Set up */ 43 /* Allocate ci array, arrays for fill computation and */ 44 /* free space for accumulating nonzero column info */ 45 ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 46 ci[0] = 0; 47 48 /* create and initialize a linked list */ 49 nlnk = bn+1; 50 ierr = PetscLLCreate(bn,bn,nlnk,lnk,lnkbt);CHKERRQ(ierr); 51 52 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ 53 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr); 54 current_space = free_space; 55 56 /* Determine symbolic info for each row of the product: */ 57 for (i=0;i<am;i++) { 58 anzi = ai[i+1] - ai[i]; 59 cnzi = 0; 60 j = anzi; 61 aj = a->j + ai[i]; 62 while (j){/* assume cols are almost in increasing order, starting from its end saves computation */ 63 j--; 64 brow = *(aj + j); 65 bnzj = bi[brow+1] - bi[brow]; 66 bjj = bj + bi[brow]; 67 /* add non-zero cols of B into the sorted linked list lnk */ 68 ierr = PetscLLAdd(bnzj,bjj,bn,nlnk,lnk,lnkbt);CHKERRQ(ierr); 69 cnzi += nlnk; 70 } 71 72 /* If free space is not available, make more free space */ 73 /* Double the amount of total space in the list */ 74 if (current_space->local_remaining<cnzi) { 75 ierr = PetscFreeSpaceGet(current_space->total_array_size,¤t_space);CHKERRQ(ierr); 76 nspacedouble++; 77 } 78 79 /* Copy data into free space, then initialize lnk */ 80 ierr = PetscLLClean(bn,bn,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 81 current_space->array += cnzi; 82 current_space->local_used += cnzi; 83 current_space->local_remaining -= cnzi; 84 85 ci[i+1] = ci[i] + cnzi; 86 } 87 88 /* Column indices are in the list of free space */ 89 /* Allocate space for cj, initialize cj, and */ 90 /* destroy list of free space and other temporary array(s) */ 91 ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 92 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 93 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 94 95 /* Allocate space for ca */ 96 ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 97 ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr); 98 99 /* put together the new symbolic matrix */ 100 ierr = MatCreateSeqAIJWithArrays(A->comm,am,bn,ci,cj,ca,C);CHKERRQ(ierr); 101 102 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 103 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 104 c = (Mat_SeqAIJ *)((*C)->data); 105 c->freedata = PETSC_TRUE; 106 c->nonew = 0; 107 108 if (nspacedouble){ 109 ierr = PetscVerboseInfo(((PetscObject)(*C),"MatMatMultSymbolic_SeqAIJ_SeqAIJ: nspacedouble:%D, nnz(A):%D, nnz(B):%D, fill:%g, nnz(C):%D\n",nspacedouble,ai[am],bi[bm],fill,ci[am]));CHKERRQ(ierr); 110 } 111 PetscFunctionReturn(0); 112 } 113 114 115 #undef __FUNCT__ 116 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqAIJ" 117 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 118 { 119 PetscErrorCode ierr; 120 PetscInt flops=0; 121 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 122 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 123 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; 124 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j; 125 PetscInt am=A->M,cm=C->M; 126 PetscInt i,j,k,anzi,bnzi,cnzi,brow,nextb; 127 MatScalar *aa=a->a,*ba=b->a,*baj,*ca=c->a; 128 129 PetscFunctionBegin; 130 /* clean old values in C */ 131 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 132 /* Traverse A row-wise. */ 133 /* Build the ith row in C by summing over nonzero columns in A, */ 134 /* the rows of B corresponding to nonzeros of A. */ 135 for (i=0;i<am;i++) { 136 anzi = ai[i+1] - ai[i]; 137 for (j=0;j<anzi;j++) { 138 brow = *aj++; 139 bnzi = bi[brow+1] - bi[brow]; 140 bjj = bj + bi[brow]; 141 baj = ba + bi[brow]; 142 nextb = 0; 143 for (k=0; nextb<bnzi; k++) { 144 if (cj[k] == bjj[nextb]){ /* ccol == bcol */ 145 ca[k] += (*aa)*baj[nextb++]; 146 } 147 } 148 flops += 2*bnzi; 149 aa++; 150 } 151 cnzi = ci[i+1] - ci[i]; 152 ca += cnzi; 153 cj += cnzi; 154 } 155 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 156 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 157 158 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatMatMultTranspose_SeqAIJ_SeqAIJ" 165 PetscErrorCode MatMatMultTranspose_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) { 166 PetscErrorCode ierr; 167 168 PetscFunctionBegin; 169 if (scall == MAT_INITIAL_MATRIX){ 170 ierr = MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 171 } 172 ierr = MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNCT__ 177 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ" 178 PetscErrorCode MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 179 { 180 PetscErrorCode ierr; 181 Mat At; 182 PetscInt *ati,*atj; 183 184 PetscFunctionBegin; 185 /* create symbolic At */ 186 ierr = MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 187 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->n,A->m,ati,atj,PETSC_NULL,&At);CHKERRQ(ierr); 188 189 /* get symbolic C=At*B */ 190 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);CHKERRQ(ierr); 191 192 /* clean up */ 193 ierr = MatDestroy(At);CHKERRQ(ierr); 194 ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 195 196 PetscFunctionReturn(0); 197 } 198 199 #undef __FUNCT__ 200 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ" 201 PetscErrorCode MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 202 { 203 PetscErrorCode ierr; 204 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data; 205 PetscInt am=A->m,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb; 206 PetscInt cm=C->m,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k,flops=0; 207 MatScalar *aa=a->a,*ba,*ca=c->a,*caj; 208 209 PetscFunctionBegin; 210 /* clear old values in C */ 211 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 212 213 /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */ 214 for (i=0;i<am;i++) { 215 bj = b->j + bi[i]; 216 ba = b->a + bi[i]; 217 bnzi = bi[i+1] - bi[i]; 218 anzi = ai[i+1] - ai[i]; 219 for (j=0; j<anzi; j++) { 220 nextb = 0; 221 crow = *aj++; 222 cjj = cj + ci[crow]; 223 caj = ca + ci[crow]; 224 /* perform sparse axpy operation. Note cjj includes bj. */ 225 for (k=0; nextb<bnzi; k++) { 226 if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */ 227 caj[k] += (*aa)*(*(ba+nextb)); 228 nextb++; 229 } 230 } 231 flops += 2*bnzi; 232 aa++; 233 } 234 } 235 236 /* Assemble the final matrix and clean up */ 237 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 238 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 239 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242