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 #define DEBUG_MATMATMULT 13 */ 14 EXTERN_C_BEGIN 15 #undef __FUNCT__ 16 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqAIJ" 17 PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 18 { 19 PetscErrorCode ierr; 20 PetscBool scalable=PETSC_FALSE,scalable_fast=PETSC_FALSE; 21 22 PetscFunctionBegin; 23 if (scall == MAT_INITIAL_MATRIX){ 24 ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr); 25 ierr = PetscOptionsBool("-matmatmult_scalable","Use a scalable but slower C=A*B","",scalable,&scalable,PETSC_NULL);CHKERRQ(ierr); 26 ierr = PetscOptionsBool("-matmatmult_scalable_fast","Use a scalable but slower C=A*B","",scalable_fast,&scalable_fast,PETSC_NULL);CHKERRQ(ierr); 27 ierr = PetscOptionsEnd();CHKERRQ(ierr); 28 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 29 if (scalable_fast){ 30 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(A,B,fill,C);CHKERRQ(ierr); 31 } else if (scalable){ 32 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,B,fill,C);CHKERRQ(ierr); 33 } else { 34 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 35 } 36 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 37 } 38 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 39 ierr = (*(*C)->ops->matmultnumeric)(A,B,*C);CHKERRQ(ierr); 40 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 EXTERN_C_END 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqAIJ" 47 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 48 { 49 PetscErrorCode ierr; 50 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 51 PetscInt *ai=a->i,*bi=b->i,*ci,*cj; 52 PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N; 53 PetscReal afill; 54 PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,nlnk_max,*lnk,ndouble=0; 55 PetscBT lnkbt; 56 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 57 58 PetscFunctionBegin; 59 /* Get ci and cj */ 60 /*---------------*/ 61 /* Allocate ci array, arrays for fill computation and */ 62 /* free space for accumulating nonzero column info */ 63 ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 64 ci[0] = 0; 65 66 /* create and initialize a linked list */ 67 nlnk_max = a->rmax*b->rmax; 68 if (!nlnk_max || nlnk_max > bn) nlnk_max = bn; 69 ierr = PetscLLCondensedCreate(nlnk_max,bn,&lnk,&lnkbt);CHKERRQ(ierr); 70 71 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ 72 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr); 73 current_space = free_space; 74 75 /* Determine ci and cj */ 76 for (i=0; i<am; i++) { 77 anzi = ai[i+1] - ai[i]; 78 aj = a->j + ai[i]; 79 for (j=0; j<anzi; j++){ 80 brow = aj[j]; 81 bnzj = bi[brow+1] - bi[brow]; 82 bj = b->j + bi[brow]; 83 /* add non-zero cols of B into the sorted linked list lnk */ 84 ierr = PetscLLCondensedAddSorted(bnzj,bj,lnk,lnkbt);CHKERRQ(ierr); 85 } 86 cnzi = lnk[0]; 87 88 /* If free space is not available, make more free space */ 89 /* Double the amount of total space in the list */ 90 if (current_space->local_remaining<cnzi) { 91 ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 92 ndouble++; 93 } 94 95 /* Copy data into free space, then initialize lnk */ 96 ierr = PetscLLCondensedClean(bn,cnzi,current_space->array,lnk,lnkbt);CHKERRQ(ierr); 97 current_space->array += cnzi; 98 current_space->local_used += cnzi; 99 current_space->local_remaining -= cnzi; 100 ci[i+1] = ci[i] + cnzi; 101 } 102 103 /* Column indices are in the list of free space */ 104 /* Allocate space for cj, initialize cj, and */ 105 /* destroy list of free space and other temporary array(s) */ 106 ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 107 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 108 ierr = PetscLLCondensedDestroy(lnk,lnkbt);CHKERRQ(ierr); 109 110 /* put together the new symbolic matrix */ 111 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bn,ci,cj,PETSC_NULL,C);CHKERRQ(ierr); 112 113 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 114 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 115 c = (Mat_SeqAIJ *)((*C)->data); 116 c->free_a = PETSC_FALSE; 117 c->free_ij = PETSC_TRUE; 118 c->nonew = 0; 119 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ; /* fast, needs non-scalable O(bn) array 'abdense' */ 120 121 /* set MatInfo */ 122 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5; 123 if (afill < 1.0) afill = 1.0; 124 c->maxnz = ci[am]; 125 c->nz = ci[am]; 126 (*C)->info.mallocs = ndouble; 127 (*C)->info.fill_ratio_given = fill; 128 (*C)->info.fill_ratio_needed = afill; 129 130 #if defined(PETSC_USE_INFO) 131 if (ci[am]) { 132 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",ndouble,fill,afill);CHKERRQ(ierr); 133 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 134 } else { 135 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 136 } 137 #endif 138 PetscFunctionReturn(0); 139 } 140 141 #undef __FUNCT__ 142 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqAIJ" 143 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 144 { 145 PetscErrorCode ierr; 146 PetscLogDouble flops=0.0; 147 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 148 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 149 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; 150 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j; 151 PetscInt am=A->rmap->n,cm=C->rmap->n; 152 PetscInt i,j,k,anzi,bnzi,cnzi,brow; 153 PetscScalar *aa=a->a,*ba=b->a,*baj,*ca,valtmp; 154 PetscScalar *ab_dense; 155 156 PetscFunctionBegin; 157 /* printf("MatMatMultNumeric_SeqAIJ_SeqAIJ...ca %p\n",c->a); */ 158 if (!c->a){ /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */ 159 ierr = PetscMalloc((ci[cm]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 160 c->a = ca; 161 c->free_a = PETSC_TRUE; 162 163 ierr = PetscMalloc(B->cmap->N*sizeof(PetscScalar),&ab_dense);CHKERRQ(ierr); 164 ierr = PetscMemzero(ab_dense,B->cmap->N*sizeof(PetscScalar));CHKERRQ(ierr); 165 c->matmult_abdense = ab_dense; 166 } else { 167 ca = c->a; 168 ab_dense = c->matmult_abdense; 169 } 170 171 /* clean old values in C */ 172 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 173 /* Traverse A row-wise. */ 174 /* Build the ith row in C by summing over nonzero columns in A, */ 175 /* the rows of B corresponding to nonzeros of A. */ 176 for (i=0; i<am; i++) { 177 anzi = ai[i+1] - ai[i]; 178 for (j=0; j<anzi; j++) { 179 brow = aj[j]; 180 bnzi = bi[brow+1] - bi[brow]; 181 bjj = bj + bi[brow]; 182 baj = ba + bi[brow]; 183 /* perform dense axpy */ 184 valtmp = aa[j]; 185 for (k=0; k<bnzi; k++) { 186 ab_dense[bjj[k]] += valtmp*baj[k]; 187 } 188 flops += 2*bnzi; 189 } 190 aj += anzi; aa += anzi; 191 192 cnzi = ci[i+1] - ci[i]; 193 for (k=0; k<cnzi; k++) { 194 ca[k] += ab_dense[cj[k]]; 195 ab_dense[cj[k]] = 0.0; /* zero ab_dense */ 196 } 197 flops += cnzi; 198 cj += cnzi; ca += cnzi; 199 } 200 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 201 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 202 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 #undef __FUNCT__ 207 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable" 208 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,Mat C) 209 { 210 PetscErrorCode ierr; 211 PetscLogDouble flops=0.0; 212 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 213 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 214 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; 215 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j; 216 PetscInt am=A->rmap->N,cm=C->rmap->N; 217 PetscInt i,j,k,anzi,bnzi,cnzi,brow; 218 PetscScalar *aa=a->a,*ba=b->a,*baj,*ca=c->a,valtmp; 219 PetscInt nextb; 220 221 PetscFunctionBegin; 222 #if defined(DEBUG_MATMATMULT) 223 ierr = PetscPrintf(PETSC_COMM_SELF,"MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable...\n");CHKERRQ(ierr); 224 #endif 225 /* clean old values in C */ 226 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 227 /* Traverse A row-wise. */ 228 /* Build the ith row in C by summing over nonzero columns in A, */ 229 /* the rows of B corresponding to nonzeros of A. */ 230 for (i=0;i<am;i++) { 231 anzi = ai[i+1] - ai[i]; 232 cnzi = ci[i+1] - ci[i]; 233 for (j=0;j<anzi;j++) { 234 brow = aj[j]; 235 bnzi = bi[brow+1] - bi[brow]; 236 bjj = bj + bi[brow]; 237 baj = ba + bi[brow]; 238 /* perform sparse axpy */ 239 valtmp = aa[j]; 240 nextb = 0; 241 for (k=0; nextb<bnzi; k++) { 242 if (cj[k] == bjj[nextb]){ /* ccol == bcol */ 243 ca[k] += valtmp*baj[nextb++]; 244 } 245 } 246 flops += 2*bnzi; 247 } 248 aj += anzi; aa += anzi; 249 cj += cnzi; ca += cnzi; 250 } 251 252 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 253 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 254 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 258 #undef __FUNCT__ 259 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast" 260 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(Mat A,Mat B,PetscReal fill,Mat *C) 261 { 262 PetscErrorCode ierr; 263 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 264 PetscInt *ai=a->i,*bi=b->i,*ci,*cj; 265 PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N; 266 MatScalar *ca; 267 PetscReal afill; 268 PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,nlnk_max,*lnk,ndouble=0; 269 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 270 271 PetscFunctionBegin; 272 /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_fast() */ 273 /*-----------------------------------------------------------------------------------------*/ 274 /* Allocate arrays for fill computation and free space for accumulating nonzero column */ 275 ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 276 ci[0] = 0; 277 278 /* create and initialize a linked list */ 279 nlnk_max = a->rmax*b->rmax; 280 if (!nlnk_max || nlnk_max > bn) nlnk_max = bn; /* in case rmax is not defined for A or B */ 281 ierr = PetscLLCondensedCreate_fast(nlnk_max,&lnk);CHKERRQ(ierr); 282 283 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ 284 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr); 285 current_space = free_space; 286 287 /* Determine ci and cj */ 288 for (i=0; i<am; i++) { 289 anzi = ai[i+1] - ai[i]; 290 aj = a->j + ai[i]; 291 for (j=0; j<anzi; j++){ 292 brow = aj[j]; 293 bnzj = bi[brow+1] - bi[brow]; 294 bj = b->j + bi[brow]; 295 /* add non-zero cols of B into the sorted linked list lnk */ 296 ierr = PetscLLCondensedAddSorted_fast(bnzj,bj,lnk);CHKERRQ(ierr); 297 } 298 cnzi = lnk[1]; 299 300 /* If free space is not available, make more free space */ 301 /* Double the amount of total space in the list */ 302 if (current_space->local_remaining<cnzi) { 303 ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 304 ndouble++; 305 } 306 307 /* Copy data into free space, then initialize lnk */ 308 ierr = PetscLLCondensedClean_fast(cnzi,current_space->array,lnk);CHKERRQ(ierr); 309 current_space->array += cnzi; 310 current_space->local_used += cnzi; 311 current_space->local_remaining -= cnzi; 312 ci[i+1] = ci[i] + cnzi; 313 } 314 315 /* Column indices are in the list of free space */ 316 /* Allocate space for cj, initialize cj, and */ 317 /* destroy list of free space and other temporary array(s) */ 318 ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 319 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 320 ierr = PetscLLCondensedDestroy_fast(lnk);CHKERRQ(ierr); 321 322 /* Allocate space for ca */ 323 ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 324 ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr); 325 326 /* put together the new symbolic matrix */ 327 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bn,ci,cj,ca,C);CHKERRQ(ierr); 328 329 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 330 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 331 c = (Mat_SeqAIJ *)((*C)->data); 332 c->free_a = PETSC_TRUE; 333 c->free_ij = PETSC_TRUE; 334 c->nonew = 0; 335 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */ 336 337 /* set MatInfo */ 338 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5; 339 if (afill < 1.0) afill = 1.0; 340 c->maxnz = ci[am]; 341 c->nz = ci[am]; 342 (*C)->info.mallocs = ndouble; 343 (*C)->info.fill_ratio_given = fill; 344 (*C)->info.fill_ratio_needed = afill; 345 346 #if defined(PETSC_USE_INFO) 347 if (ci[am]) { 348 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",ndouble,fill,afill);CHKERRQ(ierr); 349 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 350 } else { 351 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 352 } 353 #endif 354 PetscFunctionReturn(0); 355 } 356 357 358 #undef __FUNCT__ 359 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable" 360 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,PetscReal fill,Mat *C) 361 { 362 PetscErrorCode ierr; 363 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 364 PetscInt *ai=a->i,*bi=b->i,*ci,*cj; 365 PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N; 366 MatScalar *ca; 367 PetscReal afill; 368 PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,nlnk_max,*lnk,ndouble=0; 369 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 370 371 PetscFunctionBegin; 372 /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_Scalalbe() */ 373 /*---------------------------------------------------------------------------------------------*/ 374 /* Allocate arrays for fill computation and free space for accumulating nonzero column */ 375 ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 376 ci[0] = 0; 377 378 /* create and initialize a linked list */ 379 nlnk_max = a->rmax*b->rmax; 380 if (!nlnk_max || nlnk_max > bn) nlnk_max = bn; /* in case rmax is not defined for A or B */ 381 ierr = PetscLLCondensedCreate_Scalable(nlnk_max,&lnk);CHKERRQ(ierr); 382 383 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ 384 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr); 385 current_space = free_space; 386 387 /* Determine ci and cj */ 388 for (i=0; i<am; i++) { 389 anzi = ai[i+1] - ai[i]; 390 aj = a->j + ai[i]; 391 for (j=0; j<anzi; j++){ 392 brow = aj[j]; 393 bnzj = bi[brow+1] - bi[brow]; 394 bj = b->j + bi[brow]; 395 /* add non-zero cols of B into the sorted linked list lnk */ 396 ierr = PetscLLCondensedAddSorted_Scalable(bnzj,bj,lnk);CHKERRQ(ierr); 397 } 398 cnzi = lnk[0]; 399 400 /* If free space is not available, make more free space */ 401 /* Double the amount of total space in the list */ 402 if (current_space->local_remaining<cnzi) { 403 ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 404 ndouble++; 405 } 406 407 /* Copy data into free space, then initialize lnk */ 408 ierr = PetscLLCondensedClean_Scalable(cnzi,current_space->array,lnk);CHKERRQ(ierr); 409 current_space->array += cnzi; 410 current_space->local_used += cnzi; 411 current_space->local_remaining -= cnzi; 412 ci[i+1] = ci[i] + cnzi; 413 } 414 415 /* Column indices are in the list of free space */ 416 /* Allocate space for cj, initialize cj, and */ 417 /* destroy list of free space and other temporary array(s) */ 418 ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 419 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 420 ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr); 421 422 /* Allocate space for ca */ 423 /*-----------------------*/ 424 ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 425 ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr); 426 427 /* put together the new symbolic matrix */ 428 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bn,ci,cj,ca,C);CHKERRQ(ierr); 429 430 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 431 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 432 c = (Mat_SeqAIJ *)((*C)->data); 433 c->free_a = PETSC_TRUE; 434 c->free_ij = PETSC_TRUE; 435 c->nonew = 0; 436 (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */ 437 438 /* set MatInfo */ 439 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5; 440 if (afill < 1.0) afill = 1.0; 441 c->maxnz = ci[am]; 442 c->nz = ci[am]; 443 (*C)->info.mallocs = ndouble; 444 (*C)->info.fill_ratio_given = fill; 445 (*C)->info.fill_ratio_needed = afill; 446 447 #if defined(PETSC_USE_INFO) 448 if (ci[am]) { 449 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",ndouble,fill,afill);CHKERRQ(ierr); 450 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 451 } else { 452 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 453 } 454 #endif 455 PetscFunctionReturn(0); 456 } 457 458 459 /* This routine is not used. Should be removed! */ 460 #undef __FUNCT__ 461 #define __FUNCT__ "MatMatTransposeMult_SeqAIJ_SeqAIJ" 462 PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 463 { 464 PetscErrorCode ierr; 465 466 PetscFunctionBegin; 467 if (scall == MAT_INITIAL_MATRIX){ 468 ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 469 } 470 ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 471 PetscFunctionReturn(0); 472 } 473 474 #undef __FUNCT__ 475 #define __FUNCT__ "PetscContainerDestroy_Mat_MatMatTransMult" 476 PetscErrorCode PetscContainerDestroy_Mat_MatMatTransMult(void *ptr) 477 { 478 PetscErrorCode ierr; 479 Mat_MatMatTransMult *multtrans=(Mat_MatMatTransMult*)ptr; 480 481 PetscFunctionBegin; 482 ierr = MatTransposeColoringDestroy(&multtrans->matcoloring);CHKERRQ(ierr); 483 ierr = MatDestroy(&multtrans->Bt_den);CHKERRQ(ierr); 484 ierr = MatDestroy(&multtrans->ABt_den);CHKERRQ(ierr); 485 ierr = PetscFree(multtrans);CHKERRQ(ierr); 486 PetscFunctionReturn(0); 487 } 488 489 #undef __FUNCT__ 490 #define __FUNCT__ "MatDestroy_SeqAIJ_MatMatMultTrans" 491 PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(Mat A) 492 { 493 PetscErrorCode ierr; 494 PetscContainer container; 495 Mat_MatMatTransMult *multtrans=PETSC_NULL; 496 497 PetscFunctionBegin; 498 ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatTransMult",(PetscObject *)&container);CHKERRQ(ierr); 499 if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit"); 500 ierr = PetscContainerGetPointer(container,(void **)&multtrans);CHKERRQ(ierr); 501 A->ops->destroy = multtrans->destroy; 502 if (A->ops->destroy) { 503 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 504 } 505 ierr = PetscObjectCompose((PetscObject)A,"Mat_MatMatTransMult",0);CHKERRQ(ierr); 506 PetscFunctionReturn(0); 507 } 508 509 #undef __FUNCT__ 510 #define __FUNCT__ "MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ" 511 PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 512 { 513 PetscErrorCode ierr; 514 Mat Bt; 515 PetscInt *bti,*btj; 516 Mat_MatMatTransMult *multtrans; 517 PetscContainer container; 518 PetscLogDouble t0,tf,etime2=0.0; 519 520 PetscFunctionBegin; 521 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 522 /* create symbolic Bt */ 523 ierr = MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr); 524 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,B->cmap->n,B->rmap->n,bti,btj,PETSC_NULL,&Bt);CHKERRQ(ierr); 525 526 /* get symbolic C=A*Bt */ 527 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);CHKERRQ(ierr); 528 529 /* create a supporting struct for reuse intermidiate dense matrices with matcoloring */ 530 ierr = PetscNew(Mat_MatMatTransMult,&multtrans);CHKERRQ(ierr); 531 532 /* attach the supporting struct to C */ 533 ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr); 534 ierr = PetscContainerSetPointer(container,multtrans);CHKERRQ(ierr); 535 ierr = PetscContainerSetUserDestroy(container,PetscContainerDestroy_Mat_MatMatTransMult);CHKERRQ(ierr); 536 ierr = PetscObjectCompose((PetscObject)(*C),"Mat_MatMatTransMult",(PetscObject)container);CHKERRQ(ierr); 537 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 538 539 multtrans->usecoloring = PETSC_FALSE; 540 multtrans->destroy = (*C)->ops->destroy; 541 (*C)->ops->destroy = MatDestroy_SeqAIJ_MatMatMultTrans; 542 543 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 544 etime2 += tf - t0; 545 546 ierr = PetscOptionsGetBool(PETSC_NULL,"-matmattransmult_color",&multtrans->usecoloring,PETSC_NULL);CHKERRQ(ierr); 547 if (multtrans->usecoloring){ 548 /* Create MatTransposeColoring from symbolic C=A*B^T */ 549 MatTransposeColoring matcoloring; 550 ISColoring iscoloring; 551 Mat Bt_dense,C_dense; 552 PetscLogDouble etime0=0.0,etime01=0.0,etime1=0.0; 553 554 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 555 ierr = MatGetColoring(*C,MATCOLORINGLF,&iscoloring);CHKERRQ(ierr); 556 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 557 etime0 += tf - t0; 558 559 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 560 ierr = MatTransposeColoringCreate(*C,iscoloring,&matcoloring);CHKERRQ(ierr); 561 multtrans->matcoloring = matcoloring; 562 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 563 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 564 etime01 += tf - t0; 565 566 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 567 /* Create Bt_dense and C_dense = A*Bt_dense */ 568 ierr = MatCreate(PETSC_COMM_SELF,&Bt_dense);CHKERRQ(ierr); 569 ierr = MatSetSizes(Bt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);CHKERRQ(ierr); 570 ierr = MatSetType(Bt_dense,MATSEQDENSE);CHKERRQ(ierr); 571 ierr = MatSeqDenseSetPreallocation(Bt_dense,PETSC_NULL);CHKERRQ(ierr); 572 Bt_dense->assembled = PETSC_TRUE; 573 multtrans->Bt_den = Bt_dense; 574 575 ierr = MatCreate(PETSC_COMM_SELF,&C_dense);CHKERRQ(ierr); 576 ierr = MatSetSizes(C_dense,A->rmap->n,matcoloring->ncolors,A->rmap->n,matcoloring->ncolors);CHKERRQ(ierr); 577 ierr = MatSetType(C_dense,MATSEQDENSE);CHKERRQ(ierr); 578 ierr = MatSeqDenseSetPreallocation(C_dense,PETSC_NULL);CHKERRQ(ierr); 579 Bt_dense->assembled = PETSC_TRUE; 580 multtrans->ABt_den = C_dense; 581 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 582 etime1 += tf - t0; 583 584 #if defined(PETSC_USE_INFO) 585 { 586 Mat_SeqAIJ *c=(Mat_SeqAIJ*)(*C)->data; 587 ierr = PetscInfo5(*C,"Bt_dense: %D,%D; Cnz %D / (cm*ncolors %D) = %g\n",A->cmap->n,matcoloring->ncolors,c->nz,A->rmap->n*matcoloring->ncolors,(PetscReal)(c->nz)/(A->rmap->n*matcoloring->ncolors)); 588 ierr = PetscInfo5(*C,"Sym = GetColor %g + ColorCreate %g + MatDenseCreate %g + non-colorSym %g = %g\n",etime0,etime01,etime1,etime2,etime0+etime01+etime1+etime2); 589 } 590 #endif 591 } 592 /* clean up */ 593 ierr = MatDestroy(&Bt);CHKERRQ(ierr); 594 ierr = MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr); 595 596 597 598 #if defined(INEFFICIENT_ALGORITHM) 599 /* The algorithm below computes am*bm sparse inner-product - inefficient! It will be deleted later. */ 600 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 601 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 602 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*ci,*cj,*acol,*bcol; 603 PetscInt am=A->rmap->N,bm=B->rmap->N; 604 PetscInt i,j,anzi,bnzj,cnzi,nlnk,*lnk,nspacedouble=0,ka,kb,index[1]; 605 MatScalar *ca; 606 PetscBT lnkbt; 607 PetscReal afill; 608 609 /* Allocate row pointer array ci */ 610 ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 611 ci[0] = 0; 612 613 /* Create and initialize a linked list for C columns */ 614 nlnk = bm+1; 615 ierr = PetscLLCreate(bm,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr); 616 617 /* Initial FreeSpace with size fill*(nnz(A)+nnz(B)) */ 618 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr); 619 current_space = free_space; 620 621 /* Determine symbolic info for each row of the product A*B^T: */ 622 for (i=0; i<am; i++) { 623 anzi = ai[i+1] - ai[i]; 624 cnzi = 0; 625 acol = aj + ai[i]; 626 for (j=0; j<bm; j++){ 627 bnzj = bi[j+1] - bi[j]; 628 bcol= bj + bi[j]; 629 /* sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */ 630 ka = 0; kb = 0; 631 while (ka < anzi && kb < bnzj){ 632 while (acol[ka] < bcol[kb] && ka < anzi) ka++; 633 if (ka == anzi) break; 634 while (acol[ka] > bcol[kb] && kb < bnzj) kb++; 635 if (kb == bnzj) break; 636 if (acol[ka] == bcol[kb]){ /* add nonzero c(i,j) to lnk */ 637 index[0] = j; 638 ierr = PetscLLAdd(1,index,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr); 639 cnzi++; 640 break; 641 } 642 } 643 } 644 645 /* If free space is not available, make more free space */ 646 /* Double the amount of total space in the list */ 647 if (current_space->local_remaining<cnzi) { 648 ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 649 nspacedouble++; 650 } 651 652 /* Copy data into free space, then initialize lnk */ 653 ierr = PetscLLClean(bm,bm,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 654 current_space->array += cnzi; 655 current_space->local_used += cnzi; 656 current_space->local_remaining -= cnzi; 657 658 ci[i+1] = ci[i] + cnzi; 659 } 660 661 662 /* Column indices are in the list of free space. 663 Allocate array cj, copy column indices to cj, and destroy list of free space */ 664 ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 665 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 666 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 667 668 /* Allocate space for ca */ 669 ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 670 ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr); 671 672 /* put together the new symbolic matrix */ 673 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bm,ci,cj,ca,C);CHKERRQ(ierr); 674 675 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 676 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 677 c = (Mat_SeqAIJ *)((*C)->data); 678 c->free_a = PETSC_TRUE; 679 c->free_ij = PETSC_TRUE; 680 c->nonew = 0; 681 682 /* set MatInfo */ 683 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5; 684 if (afill < 1.0) afill = 1.0; 685 c->maxnz = ci[am]; 686 c->nz = ci[am]; 687 (*C)->info.mallocs = nspacedouble; 688 (*C)->info.fill_ratio_given = fill; 689 (*C)->info.fill_ratio_needed = afill; 690 691 #if defined(PETSC_USE_INFO) 692 if (ci[am]) { 693 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 694 ierr = PetscInfo1((*C),"Use MatMatTransposeMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 695 } else { 696 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 697 } 698 #endif 699 #endif 700 PetscFunctionReturn(0); 701 } 702 703 /* #define USE_ARRAY - for sparse dot product. Slower than !USE_ARRAY */ 704 #undef __FUNCT__ 705 #define __FUNCT__ "MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ" 706 PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 707 { 708 PetscErrorCode ierr; 709 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data; 710 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,anzi,bnzj,nexta,nextb,*acol,*bcol,brow; 711 PetscInt cm=C->rmap->n,*ci=c->i,*cj=c->j,i,j,cnzi,*ccol; 712 PetscLogDouble flops=0.0; 713 MatScalar *aa=a->a,*aval,*ba=b->a,*bval,*ca,*cval; 714 Mat_MatMatTransMult *multtrans; 715 PetscContainer container; 716 #if defined(USE_ARRAY) 717 MatScalar *spdot; 718 #endif 719 720 PetscFunctionBegin; 721 ierr = PetscObjectQuery((PetscObject)C,"Mat_MatMatTransMult",(PetscObject *)&container);CHKERRQ(ierr); 722 if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit"); 723 ierr = PetscContainerGetPointer(container,(void **)&multtrans);CHKERRQ(ierr); 724 if (multtrans->usecoloring){ 725 MatTransposeColoring matcoloring = multtrans->matcoloring; 726 Mat Bt_dense; 727 PetscInt m,n; 728 PetscLogDouble t0,tf,etime0=0.0,etime1=0.0,etime2=0.0; 729 Mat C_dense = multtrans->ABt_den; 730 731 Bt_dense = multtrans->Bt_den; 732 ierr = MatGetLocalSize(Bt_dense,&m,&n);CHKERRQ(ierr); 733 734 /* Get Bt_dense by Apply MatTransposeColoring to B */ 735 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 736 ierr = MatTransColoringApplySpToDen(matcoloring,B,Bt_dense);CHKERRQ(ierr); 737 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 738 etime0 += tf - t0; 739 740 /* C_dense = A*Bt_dense */ 741 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 742 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,Bt_dense,C_dense);CHKERRQ(ierr); 743 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 744 etime2 += tf - t0; 745 746 /* Recover C from C_dense */ 747 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 748 ierr = MatTransColoringApplyDenToSp(matcoloring,C_dense,C);CHKERRQ(ierr); 749 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 750 etime1 += tf - t0; 751 #if defined(PETSC_USE_INFO) 752 ierr = PetscInfo4(C,"Num = ColoringApply: %g %g + Mult_sp_dense %g = %g\n",etime0,etime1,etime2,etime0+etime1+etime2); 753 #endif 754 PetscFunctionReturn(0); 755 } 756 757 #if defined(USE_ARRAY) 758 /* allocate an array for implementing sparse inner-product */ 759 ierr = PetscMalloc((A->cmap->n+1)*sizeof(MatScalar),&spdot);CHKERRQ(ierr); 760 ierr = PetscMemzero(spdot,(A->cmap->n+1)*sizeof(MatScalar));CHKERRQ(ierr); 761 #endif 762 763 /* clear old values in C */ 764 if (!c->a){ 765 ierr = PetscMalloc((ci[cm]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 766 c->a = ca; 767 c->free_a = PETSC_TRUE; 768 } else { 769 ca = c->a; 770 } 771 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 772 773 for (i=0; i<cm; i++) { 774 anzi = ai[i+1] - ai[i]; 775 acol = aj + ai[i]; 776 aval = aa + ai[i]; 777 cnzi = ci[i+1] - ci[i]; 778 ccol = cj + ci[i]; 779 cval = ca + ci[i]; 780 for (j=0; j<cnzi; j++){ 781 brow = ccol[j]; 782 bnzj = bi[brow+1] - bi[brow]; 783 bcol = bj + bi[brow]; 784 bval = ba + bi[brow]; 785 786 /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */ 787 #if defined(USE_ARRAY) 788 /* put ba to spdot array */ 789 for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = bval[nextb]; 790 /* c(i,j)=A[i,:]*B[j,:]^T */ 791 for (nexta=0; nexta<anzi; nexta++){ 792 cval[j] += spdot[acol[nexta]]*aval[nexta]; 793 } 794 /* zero spdot array */ 795 for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = 0.0; 796 #else 797 nexta = 0; nextb = 0; 798 while (nexta<anzi && nextb<bnzj){ 799 while (acol[nexta] < bcol[nextb] && nexta < anzi) nexta++; 800 if (nexta == anzi) break; 801 while (acol[nexta] > bcol[nextb] && nextb < bnzj) nextb++; 802 if (nextb == bnzj) break; 803 if (acol[nexta] == bcol[nextb]){ 804 cval[j] += aval[nexta]*bval[nextb]; 805 nexta++; nextb++; 806 flops += 2; 807 } 808 } 809 #endif 810 } 811 } 812 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 813 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 814 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 815 #if defined(USE_ARRAY) 816 ierr = PetscFree(spdot); 817 #endif 818 PetscFunctionReturn(0); 819 } 820 821 #undef __FUNCT__ 822 #define __FUNCT__ "MatTransposeMatMult_SeqAIJ_SeqAIJ" 823 PetscErrorCode MatTransposeMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) { 824 PetscErrorCode ierr; 825 826 PetscFunctionBegin; 827 if (scall == MAT_INITIAL_MATRIX){ 828 ierr = MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 829 } 830 ierr = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 831 PetscFunctionReturn(0); 832 } 833 834 #undef __FUNCT__ 835 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ" 836 PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 837 { 838 PetscErrorCode ierr; 839 Mat At; 840 PetscInt *ati,*atj; 841 842 PetscFunctionBegin; 843 /* create symbolic At */ 844 ierr = MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 845 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->cmap->n,A->rmap->n,ati,atj,PETSC_NULL,&At);CHKERRQ(ierr); 846 847 /* get symbolic C=At*B */ 848 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);CHKERRQ(ierr); 849 850 /* clean up */ 851 ierr = MatDestroy(&At);CHKERRQ(ierr); 852 ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 853 PetscFunctionReturn(0); 854 } 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ" 858 PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 859 { 860 PetscErrorCode ierr; 861 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data; 862 PetscInt am=A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb; 863 PetscInt cm=C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k; 864 PetscLogDouble flops=0.0; 865 MatScalar *aa=a->a,*ba,*ca,*caj; 866 867 PetscFunctionBegin; 868 if (!c->a){ 869 ierr = PetscMalloc((ci[cm]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 870 c->a = ca; 871 c->free_a = PETSC_TRUE; 872 } else { 873 ca = c->a; 874 } 875 /* clear old values in C */ 876 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 877 878 /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */ 879 for (i=0;i<am;i++) { 880 bj = b->j + bi[i]; 881 ba = b->a + bi[i]; 882 bnzi = bi[i+1] - bi[i]; 883 anzi = ai[i+1] - ai[i]; 884 for (j=0; j<anzi; j++) { 885 nextb = 0; 886 crow = *aj++; 887 cjj = cj + ci[crow]; 888 caj = ca + ci[crow]; 889 /* perform sparse axpy operation. Note cjj includes bj. */ 890 for (k=0; nextb<bnzi; k++) { 891 if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */ 892 caj[k] += (*aa)*(*(ba+nextb)); 893 nextb++; 894 } 895 } 896 flops += 2*bnzi; 897 aa++; 898 } 899 } 900 901 /* Assemble the final matrix and clean up */ 902 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 903 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 904 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 905 PetscFunctionReturn(0); 906 } 907 908 EXTERN_C_BEGIN 909 #undef __FUNCT__ 910 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqDense" 911 PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 912 { 913 PetscErrorCode ierr; 914 915 PetscFunctionBegin; 916 if (scall == MAT_INITIAL_MATRIX){ 917 ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);CHKERRQ(ierr); 918 } 919 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);CHKERRQ(ierr); 920 PetscFunctionReturn(0); 921 } 922 EXTERN_C_END 923 924 #undef __FUNCT__ 925 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqDense" 926 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 927 { 928 PetscErrorCode ierr; 929 930 PetscFunctionBegin; 931 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);CHKERRQ(ierr); 932 (*C)->ops->matmult = MatMatMult_SeqAIJ_SeqDense; 933 PetscFunctionReturn(0); 934 } 935 936 #undef __FUNCT__ 937 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqDense" 938 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 939 { 940 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 941 PetscErrorCode ierr; 942 PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 943 MatScalar *aa; 944 PetscInt cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n; 945 PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam; 946 947 PetscFunctionBegin; 948 if (!cm || !cn) PetscFunctionReturn(0); 949 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); 950 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); 951 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); 952 ierr = MatGetArray(B,&b);CHKERRQ(ierr); 953 ierr = MatGetArray(C,&c);CHKERRQ(ierr); 954 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 955 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 956 colam = col*am; 957 for (i=0; i<am; i++) { /* over rows of C in those columns */ 958 r1 = r2 = r3 = r4 = 0.0; 959 n = a->i[i+1] - a->i[i]; 960 aj = a->j + a->i[i]; 961 aa = a->a + a->i[i]; 962 for (j=0; j<n; j++) { 963 r1 += (*aa)*b1[*aj]; 964 r2 += (*aa)*b2[*aj]; 965 r3 += (*aa)*b3[*aj]; 966 r4 += (*aa++)*b4[*aj++]; 967 } 968 c[colam + i] = r1; 969 c[colam + am + i] = r2; 970 c[colam + am2 + i] = r3; 971 c[colam + am3 + i] = r4; 972 } 973 b1 += bm4; 974 b2 += bm4; 975 b3 += bm4; 976 b4 += bm4; 977 } 978 for (;col<cn; col++){ /* over extra columns of C */ 979 for (i=0; i<am; i++) { /* over rows of C in those columns */ 980 r1 = 0.0; 981 n = a->i[i+1] - a->i[i]; 982 aj = a->j + a->i[i]; 983 aa = a->a + a->i[i]; 984 985 for (j=0; j<n; j++) { 986 r1 += (*aa++)*b1[*aj++]; 987 } 988 c[col*am + i] = r1; 989 } 990 b1 += bm; 991 } 992 ierr = PetscLogFlops(cn*(2.0*a->nz));CHKERRQ(ierr); 993 ierr = MatRestoreArray(B,&b);CHKERRQ(ierr); 994 ierr = MatRestoreArray(C,&c);CHKERRQ(ierr); 995 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 996 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 997 PetscFunctionReturn(0); 998 } 999 1000 /* 1001 Note very similar to MatMult_SeqAIJ(), should generate both codes from same base 1002 */ 1003 #undef __FUNCT__ 1004 #define __FUNCT__ "MatMatMultNumericAdd_SeqAIJ_SeqDense" 1005 PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 1006 { 1007 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 1008 PetscErrorCode ierr; 1009 PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 1010 MatScalar *aa; 1011 PetscInt cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n,*ii,arm; 1012 PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam,*ridx; 1013 1014 PetscFunctionBegin; 1015 if (!cm || !cn) PetscFunctionReturn(0); 1016 ierr = MatGetArray(B,&b);CHKERRQ(ierr); 1017 ierr = MatGetArray(C,&c);CHKERRQ(ierr); 1018 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 1019 1020 if (a->compressedrow.use){ /* use compressed row format */ 1021 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 1022 colam = col*am; 1023 arm = a->compressedrow.nrows; 1024 ii = a->compressedrow.i; 1025 ridx = a->compressedrow.rindex; 1026 for (i=0; i<arm; i++) { /* over rows of C in those columns */ 1027 r1 = r2 = r3 = r4 = 0.0; 1028 n = ii[i+1] - ii[i]; 1029 aj = a->j + ii[i]; 1030 aa = a->a + ii[i]; 1031 for (j=0; j<n; j++) { 1032 r1 += (*aa)*b1[*aj]; 1033 r2 += (*aa)*b2[*aj]; 1034 r3 += (*aa)*b3[*aj]; 1035 r4 += (*aa++)*b4[*aj++]; 1036 } 1037 c[colam + ridx[i]] += r1; 1038 c[colam + am + ridx[i]] += r2; 1039 c[colam + am2 + ridx[i]] += r3; 1040 c[colam + am3 + ridx[i]] += r4; 1041 } 1042 b1 += bm4; 1043 b2 += bm4; 1044 b3 += bm4; 1045 b4 += bm4; 1046 } 1047 for (;col<cn; col++){ /* over extra columns of C */ 1048 colam = col*am; 1049 arm = a->compressedrow.nrows; 1050 ii = a->compressedrow.i; 1051 ridx = a->compressedrow.rindex; 1052 for (i=0; i<arm; i++) { /* over rows of C in those columns */ 1053 r1 = 0.0; 1054 n = ii[i+1] - ii[i]; 1055 aj = a->j + ii[i]; 1056 aa = a->a + ii[i]; 1057 1058 for (j=0; j<n; j++) { 1059 r1 += (*aa++)*b1[*aj++]; 1060 } 1061 c[col*am + ridx[i]] += r1; 1062 } 1063 b1 += bm; 1064 } 1065 } else { 1066 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 1067 colam = col*am; 1068 for (i=0; i<am; i++) { /* over rows of C in those columns */ 1069 r1 = r2 = r3 = r4 = 0.0; 1070 n = a->i[i+1] - a->i[i]; 1071 aj = a->j + a->i[i]; 1072 aa = a->a + a->i[i]; 1073 for (j=0; j<n; j++) { 1074 r1 += (*aa)*b1[*aj]; 1075 r2 += (*aa)*b2[*aj]; 1076 r3 += (*aa)*b3[*aj]; 1077 r4 += (*aa++)*b4[*aj++]; 1078 } 1079 c[colam + i] += r1; 1080 c[colam + am + i] += r2; 1081 c[colam + am2 + i] += r3; 1082 c[colam + am3 + i] += r4; 1083 } 1084 b1 += bm4; 1085 b2 += bm4; 1086 b3 += bm4; 1087 b4 += bm4; 1088 } 1089 for (;col<cn; col++){ /* over extra columns of C */ 1090 for (i=0; i<am; i++) { /* over rows of C in those columns */ 1091 r1 = 0.0; 1092 n = a->i[i+1] - a->i[i]; 1093 aj = a->j + a->i[i]; 1094 aa = a->a + a->i[i]; 1095 1096 for (j=0; j<n; j++) { 1097 r1 += (*aa++)*b1[*aj++]; 1098 } 1099 c[col*am + i] += r1; 1100 } 1101 b1 += bm; 1102 } 1103 } 1104 ierr = PetscLogFlops(cn*2.0*a->nz);CHKERRQ(ierr); 1105 ierr = MatRestoreArray(B,&b);CHKERRQ(ierr); 1106 ierr = MatRestoreArray(C,&c);CHKERRQ(ierr); 1107 PetscFunctionReturn(0); 1108 } 1109 1110 #undef __FUNCT__ 1111 #define __FUNCT__ "MatTransColoringApplySpToDen_SeqAIJ" 1112 PetscErrorCode MatTransColoringApplySpToDen_SeqAIJ(MatTransposeColoring coloring,Mat B,Mat Btdense) 1113 { 1114 PetscErrorCode ierr; 1115 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1116 Mat_SeqDense *btdense = (Mat_SeqDense*)Btdense->data; 1117 PetscInt *bi=b->i,*bj=b->j; 1118 PetscInt m=Btdense->rmap->n,n=Btdense->cmap->n,j,k,l,col,anz,*btcol,brow,ncolumns; 1119 MatScalar *btval,*btval_den,*ba=b->a; 1120 PetscInt *columns=coloring->columns,*colorforcol=coloring->colorforcol,ncolors=coloring->ncolors; 1121 1122 PetscFunctionBegin; 1123 btval_den=btdense->v; 1124 ierr = PetscMemzero(btval_den,(m*n)*sizeof(MatScalar));CHKERRQ(ierr); 1125 for (k=0; k<ncolors; k++) { 1126 ncolumns = coloring->ncolumns[k]; 1127 for (l=0; l<ncolumns; l++) { /* insert a row of B to a column of Btdense */ 1128 col = *(columns + colorforcol[k] + l); 1129 btcol = bj + bi[col]; 1130 btval = ba + bi[col]; 1131 anz = bi[col+1] - bi[col]; 1132 for (j=0; j<anz; j++){ 1133 brow = btcol[j]; 1134 btval_den[brow] = btval[j]; 1135 } 1136 } 1137 btval_den += m; 1138 } 1139 PetscFunctionReturn(0); 1140 } 1141 1142 #undef __FUNCT__ 1143 #define __FUNCT__ "MatTransColoringApplyDenToSp_SeqAIJ" 1144 PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring,Mat Cden,Mat Csp) 1145 { 1146 PetscErrorCode ierr; 1147 Mat_SeqAIJ *csp = (Mat_SeqAIJ*)Csp->data; 1148 PetscInt k,l,*row,*idx,m,ncolors=matcoloring->ncolors,nrows; 1149 PetscScalar *ca_den,*cp_den,*ca=csp->a; 1150 PetscInt *rows=matcoloring->rows,*spidx=matcoloring->columnsforspidx,*colorforrow=matcoloring->colorforrow; 1151 1152 PetscFunctionBegin; 1153 ierr = MatGetLocalSize(Csp,&m,PETSC_NULL);CHKERRQ(ierr); 1154 ierr = MatGetArray(Cden,&ca_den);CHKERRQ(ierr); 1155 cp_den = ca_den; 1156 for (k=0; k<ncolors; k++) { 1157 nrows = matcoloring->nrows[k]; 1158 row = rows + colorforrow[k]; 1159 idx = spidx + colorforrow[k]; 1160 for (l=0; l<nrows; l++){ 1161 ca[idx[l]] = cp_den[row[l]]; 1162 } 1163 cp_den += m; 1164 } 1165 ierr = MatRestoreArray(Cden,&ca_den);CHKERRQ(ierr); 1166 PetscFunctionReturn(0); 1167 } 1168 1169 /* 1170 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 1171 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 1172 spidx[], index of a->j, to be used for setting 'columnsforspidx' in MatTransposeColoringCreate_SeqAIJ(). 1173 */ 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ_Color" 1176 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 1177 { 1178 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1179 PetscErrorCode ierr; 1180 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 1181 PetscInt nz = a->i[m],row,*jj,mr,col; 1182 PetscInt *cspidx; 1183 1184 PetscFunctionBegin; 1185 *nn = n; 1186 if (!ia) PetscFunctionReturn(0); 1187 if (symmetric) { 1188 SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"MatGetColumnIJ_SeqAIJ_Color() not supported for the case symmetric"); 1189 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 1190 } else { 1191 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 1192 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 1193 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 1194 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 1195 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cspidx);CHKERRQ(ierr); 1196 jj = a->j; 1197 for (i=0; i<nz; i++) { 1198 collengths[jj[i]]++; 1199 } 1200 cia[0] = oshift; 1201 for (i=0; i<n; i++) { 1202 cia[i+1] = cia[i] + collengths[i]; 1203 } 1204 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 1205 jj = a->j; 1206 for (row=0; row<m; row++) { 1207 mr = a->i[row+1] - a->i[row]; 1208 for (i=0; i<mr; i++) { 1209 col = *jj++; 1210 cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */ 1211 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 1212 } 1213 } 1214 ierr = PetscFree(collengths);CHKERRQ(ierr); 1215 *ia = cia; *ja = cja; 1216 *spidx = cspidx; 1217 } 1218 PetscFunctionReturn(0); 1219 } 1220 1221 #undef __FUNCT__ 1222 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ_Color" 1223 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 1224 { 1225 PetscErrorCode ierr; 1226 1227 PetscFunctionBegin; 1228 if (!ia) PetscFunctionReturn(0); 1229 1230 ierr = PetscFree(*ia);CHKERRQ(ierr); 1231 ierr = PetscFree(*ja);CHKERRQ(ierr); 1232 ierr = PetscFree(*spidx);CHKERRQ(ierr); 1233 PetscFunctionReturn(0); 1234 } 1235 1236 #undef __FUNCT__ 1237 #define __FUNCT__ "MatTransposeColoringCreate_SeqAIJ" 1238 PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatTransposeColoring c) 1239 { 1240 PetscErrorCode ierr; 1241 PetscInt i,n,nrows,N,j,k,m,*row_idx,*ci,*cj,ncols,col,cm; 1242 const PetscInt *is; 1243 PetscInt nis = iscoloring->n,*rowhit,bs = 1; 1244 IS *isa; 1245 PetscBool done; 1246 PetscBool flg1,flg2; 1247 Mat_SeqAIJ *csp = (Mat_SeqAIJ*)mat->data; 1248 PetscInt *colorforrow,*rows,*rows_i,*columnsforspidx,*columnsforspidx_i,*idxhit,*spidx; 1249 PetscInt *colorforcol,*columns,*columns_i; 1250 1251 PetscFunctionBegin; 1252 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 1253 1254 /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */ 1255 ierr = PetscTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);CHKERRQ(ierr); 1256 ierr = PetscTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);CHKERRQ(ierr); 1257 if (flg1 || flg2) { 1258 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1259 } 1260 1261 N = mat->cmap->N/bs; 1262 c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */ 1263 c->N = mat->cmap->N/bs; 1264 c->m = mat->rmap->N/bs; 1265 c->rstart = 0; 1266 1267 c->ncolors = nis; 1268 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 1269 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 1270 ierr = PetscMalloc2(csp->nz+1,PetscInt,&rows,csp->nz+1,PetscInt,&columnsforspidx);CHKERRQ(ierr); 1271 ierr = PetscMalloc((nis+1)*sizeof(PetscInt),&colorforrow);CHKERRQ(ierr); 1272 colorforrow[0] = 0; 1273 rows_i = rows; 1274 columnsforspidx_i = columnsforspidx; 1275 1276 ierr = PetscMalloc((nis+1)*sizeof(PetscInt),&colorforcol);CHKERRQ(ierr); 1277 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&columns);CHKERRQ(ierr); 1278 colorforcol[0] = 0; 1279 columns_i = columns; 1280 1281 ierr = MatGetColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,&done);CHKERRQ(ierr); /* column-wise storage of mat */ 1282 if (!done) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"MatGetColumnIJ() not supported for matrix type %s",((PetscObject)mat)->type_name); 1283 1284 cm = c->m; 1285 ierr = PetscMalloc((cm+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 1286 ierr = PetscMalloc((cm+1)*sizeof(PetscInt),&idxhit);CHKERRQ(ierr); 1287 for (i=0; i<nis; i++) { 1288 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 1289 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 1290 c->ncolumns[i] = n; 1291 if (n) { 1292 ierr = PetscMemcpy(columns_i,is,n*sizeof(PetscInt));CHKERRQ(ierr); 1293 } 1294 colorforcol[i+1] = colorforcol[i] + n; 1295 columns_i += n; 1296 1297 /* fast, crude version requires O(N*N) work */ 1298 ierr = PetscMemzero(rowhit,cm*sizeof(PetscInt));CHKERRQ(ierr); 1299 1300 /* loop over columns*/ 1301 for (j=0; j<n; j++) { 1302 col = is[j]; 1303 row_idx = cj + ci[col]; 1304 m = ci[col+1] - ci[col]; 1305 /* loop over columns marking them in rowhit */ 1306 for (k=0; k<m; k++) { 1307 idxhit[*row_idx] = spidx[ci[col] + k]; 1308 rowhit[*row_idx++] = col + 1; 1309 } 1310 } 1311 /* count the number of hits */ 1312 nrows = 0; 1313 for (j=0; j<cm; j++) { 1314 if (rowhit[j]) nrows++; 1315 } 1316 c->nrows[i] = nrows; 1317 colorforrow[i+1] = colorforrow[i] + nrows; 1318 1319 nrows = 0; 1320 for (j=0; j<cm; j++) { 1321 if (rowhit[j]) { 1322 rows_i[nrows] = j; 1323 columnsforspidx_i[nrows] = idxhit[j]; 1324 nrows++; 1325 } 1326 } 1327 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 1328 rows_i += nrows; columnsforspidx_i += nrows; 1329 } 1330 ierr = MatRestoreColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,&done);CHKERRQ(ierr); 1331 ierr = PetscFree(rowhit);CHKERRQ(ierr); 1332 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 1333 #if defined(PETSC_USE_DEBUG) 1334 if (csp->nz != colorforrow[nis]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"csp->nz %d != colorforrow[nis] %d",csp->nz,colorforrow[nis]); 1335 #endif 1336 1337 c->colorforrow = colorforrow; 1338 c->rows = rows; 1339 c->columnsforspidx = columnsforspidx; 1340 c->colorforcol = colorforcol; 1341 c->columns = columns; 1342 ierr = PetscFree(idxhit);CHKERRQ(ierr); 1343 PetscFunctionReturn(0); 1344 } 1345