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 <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 11 12 EXTERN_C_BEGIN 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 EXTERN_C_END 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->rmap->N,bn=B->cmap->N,bm=B->rmap->N; 37 PetscInt i,j,anzi,brow,bnzj,cnzi,nlnk,*lnk,nspacedouble=0; 38 MatScalar *ca; 39 PetscBT lnkbt; 40 PetscReal afill; 41 42 PetscFunctionBegin; 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(cnzi+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(((PetscObject)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->free_a = PETSC_TRUE; 106 c->free_ij = PETSC_TRUE; 107 c->nonew = 0; 108 109 /* set MatInfo */ 110 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5; 111 if (afill < 1.0) afill = 1.0; 112 c->maxnz = ci[am]; 113 c->nz = ci[am]; 114 (*C)->info.mallocs = nspacedouble; 115 (*C)->info.fill_ratio_given = fill; 116 (*C)->info.fill_ratio_needed = afill; 117 118 #if defined(PETSC_USE_INFO) 119 if (ci[am]) { 120 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 121 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 122 } else { 123 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 124 } 125 #endif 126 PetscFunctionReturn(0); 127 } 128 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqAIJ" 132 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 133 { 134 PetscErrorCode ierr; 135 PetscLogDouble flops=0.0; 136 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 137 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 138 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; 139 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j; 140 PetscInt am=A->rmap->N,cm=C->rmap->N; 141 PetscInt i,j,k,anzi,bnzi,cnzi,brow,nextb; 142 MatScalar *aa=a->a,*ba=b->a,*baj,*ca=c->a; 143 144 PetscFunctionBegin; 145 /* clean old values in C */ 146 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 147 /* Traverse A row-wise. */ 148 /* Build the ith row in C by summing over nonzero columns in A, */ 149 /* the rows of B corresponding to nonzeros of A. */ 150 for (i=0;i<am;i++) { 151 anzi = ai[i+1] - ai[i]; 152 for (j=0;j<anzi;j++) { 153 brow = *aj++; 154 bnzi = bi[brow+1] - bi[brow]; 155 bjj = bj + bi[brow]; 156 baj = ba + bi[brow]; 157 nextb = 0; 158 for (k=0; nextb<bnzi; k++) { 159 if (cj[k] == bjj[nextb]){ /* ccol == bcol */ 160 ca[k] += (*aa)*baj[nextb++]; 161 } 162 } 163 flops += 2*bnzi; 164 aa++; 165 } 166 cnzi = ci[i+1] - ci[i]; 167 ca += cnzi; 168 cj += cnzi; 169 } 170 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 171 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 172 173 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 174 PetscFunctionReturn(0); 175 } 176 177 #undef __FUNCT__ 178 #define __FUNCT__ "MatMatMultTranspose_SeqAIJ_SeqAIJ" 179 PetscErrorCode MatMatMultTranspose_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 180 { 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 if (scall == MAT_INITIAL_MATRIX){ 185 ierr = MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 186 } 187 ierr = MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 #undef __FUNCT__ 192 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ" 193 PetscErrorCode MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 194 { 195 PetscErrorCode ierr; 196 Mat Bt; 197 PetscInt *bti,*btj; 198 199 PetscFunctionBegin; 200 /* create symbolic Bt */ 201 ierr = MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr); 202 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,B->cmap->n,B->rmap->n,bti,btj,PETSC_NULL,&Bt);CHKERRQ(ierr); 203 204 /* get symbolic C=A*Bt */ 205 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);CHKERRQ(ierr); 206 207 /* clean up */ 208 ierr = MatDestroy(&Bt);CHKERRQ(ierr); 209 ierr = MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr); 210 211 #if defined(INEFFICIENT_ALGORITHM) 212 /* The algorithm below computes am*bm sparse inner-product - inefficient! It will be deleted later. */ 213 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 214 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 215 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*ci,*cj,*acol,*bcol; 216 PetscInt am=A->rmap->N,bm=B->rmap->N; 217 PetscInt i,j,anzi,bnzj,cnzi,nlnk,*lnk,nspacedouble=0,ka,kb,index[1]; 218 MatScalar *ca; 219 PetscBT lnkbt; 220 PetscReal afill; 221 222 /* Allocate row pointer array ci */ 223 ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 224 ci[0] = 0; 225 226 /* Create and initialize a linked list for C columns */ 227 nlnk = bm+1; 228 ierr = PetscLLCreate(bm,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr); 229 230 /* Initial FreeSpace with size fill*(nnz(A)+nnz(B)) */ 231 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr); 232 current_space = free_space; 233 234 /* Determine symbolic info for each row of the product A*B^T: */ 235 for (i=0; i<am; i++) { 236 anzi = ai[i+1] - ai[i]; 237 cnzi = 0; 238 acol = aj + ai[i]; 239 for (j=0; j<bm; j++){ 240 bnzj = bi[j+1] - bi[j]; 241 bcol= bj + bi[j]; 242 /* sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */ 243 ka = 0; kb = 0; 244 while (ka < anzi && kb < bnzj){ 245 while (acol[ka] < bcol[kb] && ka < anzi) ka++; 246 if (ka == anzi) break; 247 while (acol[ka] > bcol[kb] && kb < bnzj) kb++; 248 if (kb == bnzj) break; 249 if (acol[ka] == bcol[kb]){ /* add nonzero c(i,j) to lnk */ 250 index[0] = j; 251 ierr = PetscLLAdd(1,index,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr); 252 cnzi++; 253 break; 254 } 255 } 256 } 257 258 /* If free space is not available, make more free space */ 259 /* Double the amount of total space in the list */ 260 if (current_space->local_remaining<cnzi) { 261 ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 262 nspacedouble++; 263 } 264 265 /* Copy data into free space, then initialize lnk */ 266 ierr = PetscLLClean(bm,bm,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 267 current_space->array += cnzi; 268 current_space->local_used += cnzi; 269 current_space->local_remaining -= cnzi; 270 271 ci[i+1] = ci[i] + cnzi; 272 } 273 274 275 /* Column indices are in the list of free space. 276 Allocate array cj, copy column indices to cj, and destroy list of free space */ 277 ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 278 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 279 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 280 281 /* Allocate space for ca */ 282 ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 283 ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr); 284 285 /* put together the new symbolic matrix */ 286 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bm,ci,cj,ca,C);CHKERRQ(ierr); 287 288 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 289 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 290 c = (Mat_SeqAIJ *)((*C)->data); 291 c->free_a = PETSC_TRUE; 292 c->free_ij = PETSC_TRUE; 293 c->nonew = 0; 294 295 /* set MatInfo */ 296 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5; 297 if (afill < 1.0) afill = 1.0; 298 c->maxnz = ci[am]; 299 c->nz = ci[am]; 300 (*C)->info.mallocs = nspacedouble; 301 (*C)->info.fill_ratio_given = fill; 302 (*C)->info.fill_ratio_needed = afill; 303 304 #if defined(PETSC_USE_INFO) 305 if (ci[am]) { 306 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 307 ierr = PetscInfo1((*C),"Use MatMatMultTranspose(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 308 } else { 309 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 310 } 311 #endif 312 #endif 313 PetscFunctionReturn(0); 314 } 315 316 /* #define USE_ARRAY - for sparse dot product. Slower than !USE_ARRAY */ 317 #undef __FUNCT__ 318 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ" 319 PetscErrorCode MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 320 { 321 PetscErrorCode ierr; 322 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data; 323 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,anzi,bnzj,nexta,nextb,*acol,*bcol,brow; 324 PetscInt cm=C->rmap->n,*ci=c->i,*cj=c->j,i,j,cnzi,*ccol; 325 PetscLogDouble flops=0.0; 326 MatScalar *aa=a->a,*aval,*ba=b->a,*bval,*ca=c->a,*cval; 327 #if defined(USE_ARRAY) 328 MatScalar *spdot; 329 #endif 330 331 PetscFunctionBegin; 332 #if defined(USE_ARRAY) 333 /* allocate an array for implementing sparse inner-product */ 334 ierr = PetscMalloc((A->cmap->n+1)*sizeof(MatScalar),&spdot);CHKERRQ(ierr); 335 ierr = PetscMemzero(spdot,(A->cmap->n+1)*sizeof(MatScalar));CHKERRQ(ierr); 336 #endif 337 338 /* clear old values in C */ 339 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 340 341 for (i=0; i<cm; i++) { 342 anzi = ai[i+1] - ai[i]; 343 acol = aj + ai[i]; 344 aval = aa + ai[i]; 345 cnzi = ci[i+1] - ci[i]; 346 ccol = cj + ci[i]; 347 cval = ca + ci[i]; 348 for (j=0; j<cnzi; j++){ 349 brow = ccol[j]; 350 bnzj = bi[brow+1] - bi[brow]; 351 bcol = bj + bi[brow]; 352 bval = ba + bi[brow]; 353 354 /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */ 355 #if defined(USE_ARRAY) 356 /* put ba to spdot array */ 357 for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = bval[nextb]; 358 /* c(i,j)=A[i,:]*B[j,:]^T */ 359 for (nexta=0; nexta<anzi; nexta++){ 360 cval[j] += spdot[acol[nexta]]*aval[nexta]; 361 } 362 /* zero spdot array */ 363 for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = 0.0; 364 #else 365 nexta = 0; nextb = 0; 366 while (nexta<anzi && nextb<bnzj){ 367 while (acol[nexta] < bcol[nextb] && nexta < anzi) nexta++; 368 if (nexta == anzi) break; 369 while (acol[nexta] > bcol[nextb] && nextb < bnzj) nextb++; 370 if (nextb == bnzj) break; 371 if (acol[nexta] == bcol[nextb]){ 372 cval[j] += aval[nexta]*bval[nextb]; 373 nexta++; nextb++; 374 flops += 2; 375 } 376 } 377 #endif 378 } 379 } 380 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 381 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 382 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 383 #if defined(USE_ARRAY) 384 ierr = PetscFree(spdot); 385 #endif 386 PetscFunctionReturn(0); 387 } 388 389 #undef __FUNCT__ 390 #define __FUNCT__ "MatMatTransposeMult_SeqAIJ_SeqAIJ" 391 PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) { 392 PetscErrorCode ierr; 393 394 PetscFunctionBegin; 395 if (scall == MAT_INITIAL_MATRIX){ 396 ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 397 } 398 ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 399 PetscFunctionReturn(0); 400 } 401 402 #undef __FUNCT__ 403 #define __FUNCT__ "MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ" 404 PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 405 { 406 PetscErrorCode ierr; 407 Mat At; 408 PetscInt *ati,*atj; 409 410 PetscFunctionBegin; 411 /* create symbolic At */ 412 ierr = MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 413 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->cmap->n,A->rmap->n,ati,atj,PETSC_NULL,&At);CHKERRQ(ierr); 414 415 /* get symbolic C=At*B */ 416 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);CHKERRQ(ierr); 417 418 /* clean up */ 419 ierr = MatDestroy(&At);CHKERRQ(ierr); 420 ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 421 PetscFunctionReturn(0); 422 } 423 424 #undef __FUNCT__ 425 #define __FUNCT__ "MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ" 426 PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 427 { 428 PetscErrorCode ierr; 429 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data; 430 PetscInt am=A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb; 431 PetscInt cm=C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k; 432 PetscLogDouble flops=0.0; 433 MatScalar *aa=a->a,*ba,*ca=c->a,*caj; 434 435 PetscFunctionBegin; 436 /* clear old values in C */ 437 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 438 439 /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */ 440 for (i=0;i<am;i++) { 441 bj = b->j + bi[i]; 442 ba = b->a + bi[i]; 443 bnzi = bi[i+1] - bi[i]; 444 anzi = ai[i+1] - ai[i]; 445 for (j=0; j<anzi; j++) { 446 nextb = 0; 447 crow = *aj++; 448 cjj = cj + ci[crow]; 449 caj = ca + ci[crow]; 450 /* perform sparse axpy operation. Note cjj includes bj. */ 451 for (k=0; nextb<bnzi; k++) { 452 if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */ 453 caj[k] += (*aa)*(*(ba+nextb)); 454 nextb++; 455 } 456 } 457 flops += 2*bnzi; 458 aa++; 459 } 460 } 461 462 /* Assemble the final matrix and clean up */ 463 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 464 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 465 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 466 PetscFunctionReturn(0); 467 } 468 469 EXTERN_C_BEGIN 470 #undef __FUNCT__ 471 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqDense" 472 PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 473 { 474 PetscErrorCode ierr; 475 476 PetscFunctionBegin; 477 if (scall == MAT_INITIAL_MATRIX){ 478 ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);CHKERRQ(ierr); 479 } 480 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);CHKERRQ(ierr); 481 PetscFunctionReturn(0); 482 } 483 EXTERN_C_END 484 485 #undef __FUNCT__ 486 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqDense" 487 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 488 { 489 PetscErrorCode ierr; 490 491 PetscFunctionBegin; 492 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);CHKERRQ(ierr); 493 PetscFunctionReturn(0); 494 } 495 496 #undef __FUNCT__ 497 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqDense" 498 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 499 { 500 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 501 PetscErrorCode ierr; 502 PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 503 MatScalar *aa; 504 PetscInt cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n; 505 PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam; 506 507 PetscFunctionBegin; 508 if (!cm || !cn) PetscFunctionReturn(0); 509 if (bm != 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,bm); 510 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); 511 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); 512 ierr = MatGetArray(B,&b);CHKERRQ(ierr); 513 ierr = MatGetArray(C,&c);CHKERRQ(ierr); 514 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 515 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 516 colam = col*am; 517 for (i=0; i<am; i++) { /* over rows of C in those columns */ 518 r1 = r2 = r3 = r4 = 0.0; 519 n = a->i[i+1] - a->i[i]; 520 aj = a->j + a->i[i]; 521 aa = a->a + a->i[i]; 522 for (j=0; j<n; j++) { 523 r1 += (*aa)*b1[*aj]; 524 r2 += (*aa)*b2[*aj]; 525 r3 += (*aa)*b3[*aj]; 526 r4 += (*aa++)*b4[*aj++]; 527 } 528 c[colam + i] = r1; 529 c[colam + am + i] = r2; 530 c[colam + am2 + i] = r3; 531 c[colam + am3 + i] = r4; 532 } 533 b1 += bm4; 534 b2 += bm4; 535 b3 += bm4; 536 b4 += bm4; 537 } 538 for (;col<cn; col++){ /* over extra columns of C */ 539 for (i=0; i<am; i++) { /* over rows of C in those columns */ 540 r1 = 0.0; 541 n = a->i[i+1] - a->i[i]; 542 aj = a->j + a->i[i]; 543 aa = a->a + a->i[i]; 544 545 for (j=0; j<n; j++) { 546 r1 += (*aa++)*b1[*aj++]; 547 } 548 c[col*am + i] = r1; 549 } 550 b1 += bm; 551 } 552 ierr = PetscLogFlops(cn*(2.0*a->nz));CHKERRQ(ierr); 553 ierr = MatRestoreArray(B,&b);CHKERRQ(ierr); 554 ierr = MatRestoreArray(C,&c);CHKERRQ(ierr); 555 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 556 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 557 PetscFunctionReturn(0); 558 } 559 560 /* 561 Note very similar to MatMult_SeqAIJ(), should generate both codes from same base 562 */ 563 #undef __FUNCT__ 564 #define __FUNCT__ "MatMatMultNumericAdd_SeqAIJ_SeqDense" 565 PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 566 { 567 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 568 PetscErrorCode ierr; 569 PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 570 MatScalar *aa; 571 PetscInt cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n,*ii,arm; 572 PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam,*ridx; 573 574 PetscFunctionBegin; 575 if (!cm || !cn) PetscFunctionReturn(0); 576 ierr = MatGetArray(B,&b);CHKERRQ(ierr); 577 ierr = MatGetArray(C,&c);CHKERRQ(ierr); 578 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 579 580 if (a->compressedrow.use){ /* use compressed row format */ 581 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 582 colam = col*am; 583 arm = a->compressedrow.nrows; 584 ii = a->compressedrow.i; 585 ridx = a->compressedrow.rindex; 586 for (i=0; i<arm; i++) { /* over rows of C in those columns */ 587 r1 = r2 = r3 = r4 = 0.0; 588 n = ii[i+1] - ii[i]; 589 aj = a->j + ii[i]; 590 aa = a->a + ii[i]; 591 for (j=0; j<n; j++) { 592 r1 += (*aa)*b1[*aj]; 593 r2 += (*aa)*b2[*aj]; 594 r3 += (*aa)*b3[*aj]; 595 r4 += (*aa++)*b4[*aj++]; 596 } 597 c[colam + ridx[i]] += r1; 598 c[colam + am + ridx[i]] += r2; 599 c[colam + am2 + ridx[i]] += r3; 600 c[colam + am3 + ridx[i]] += r4; 601 } 602 b1 += bm4; 603 b2 += bm4; 604 b3 += bm4; 605 b4 += bm4; 606 } 607 for (;col<cn; col++){ /* over extra columns of C */ 608 colam = col*am; 609 arm = a->compressedrow.nrows; 610 ii = a->compressedrow.i; 611 ridx = a->compressedrow.rindex; 612 for (i=0; i<arm; i++) { /* over rows of C in those columns */ 613 r1 = 0.0; 614 n = ii[i+1] - ii[i]; 615 aj = a->j + ii[i]; 616 aa = a->a + ii[i]; 617 618 for (j=0; j<n; j++) { 619 r1 += (*aa++)*b1[*aj++]; 620 } 621 c[col*am + ridx[i]] += r1; 622 } 623 b1 += bm; 624 } 625 } else { 626 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 627 colam = col*am; 628 for (i=0; i<am; i++) { /* over rows of C in those columns */ 629 r1 = r2 = r3 = r4 = 0.0; 630 n = a->i[i+1] - a->i[i]; 631 aj = a->j + a->i[i]; 632 aa = a->a + a->i[i]; 633 for (j=0; j<n; j++) { 634 r1 += (*aa)*b1[*aj]; 635 r2 += (*aa)*b2[*aj]; 636 r3 += (*aa)*b3[*aj]; 637 r4 += (*aa++)*b4[*aj++]; 638 } 639 c[colam + i] += r1; 640 c[colam + am + i] += r2; 641 c[colam + am2 + i] += r3; 642 c[colam + am3 + i] += r4; 643 } 644 b1 += bm4; 645 b2 += bm4; 646 b3 += bm4; 647 b4 += bm4; 648 } 649 for (;col<cn; col++){ /* over extra columns of C */ 650 for (i=0; i<am; i++) { /* over rows of C in those columns */ 651 r1 = 0.0; 652 n = a->i[i+1] - a->i[i]; 653 aj = a->j + a->i[i]; 654 aa = a->a + a->i[i]; 655 656 for (j=0; j<n; j++) { 657 r1 += (*aa++)*b1[*aj++]; 658 } 659 c[col*am + i] += r1; 660 } 661 b1 += bm; 662 } 663 } 664 ierr = PetscLogFlops(cn*2.0*a->nz);CHKERRQ(ierr); 665 ierr = MatRestoreArray(B,&b);CHKERRQ(ierr); 666 ierr = MatRestoreArray(C,&c);CHKERRQ(ierr); 667 PetscFunctionReturn(0); 668 } 669