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