1 2 /* 3 Defines matrix-matrix product routines for pairs of SeqAIJ matrices 4 C = A * B 5 */ 6 7 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 8 #include <../src/mat/utils/freespace.h> 9 #include <petscbt.h> 10 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 11 12 EXTERN_C_BEGIN 13 #undef __FUNCT__ 14 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqAIJ" 15 PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 16 { 17 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 if (scall == MAT_INITIAL_MATRIX){ 21 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 22 } 23 ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 24 PetscFunctionReturn(0); 25 } 26 EXTERN_C_END 27 28 #undef __FUNCT__ 29 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqAIJ" 30 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 31 { 32 PetscErrorCode ierr; 33 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 34 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 35 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci,*cj; 36 PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N; 37 PetscInt i,j,anzi,brow,bnzj,cnzi,nlnk,*lnk,nspacedouble=0; 38 MatScalar *ca; 39 PetscBT lnkbt; 40 PetscReal afill; 41 42 PetscFunctionBegin; 43 /* Allocate ci array, arrays for fill computation and */ 44 /* free space for accumulating nonzero column info */ 45 ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 46 ci[0] = 0; 47 48 /* create and initialize a linked list */ 49 nlnk = bn+1; 50 ierr = PetscLLCreate(bn,bn,nlnk,lnk,lnkbt);CHKERRQ(ierr); 51 52 /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ 53 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr); 54 current_space = free_space; 55 56 /* Determine symbolic info for each row of the product: */ 57 for (i=0;i<am;i++) { 58 anzi = ai[i+1] - ai[i]; 59 cnzi = 0; 60 j = anzi; 61 aj = a->j + ai[i]; 62 while (j){/* assume cols are almost in increasing order, starting from its end saves computation */ 63 j--; 64 brow = *(aj + j); 65 bnzj = bi[brow+1] - bi[brow]; 66 bjj = bj + bi[brow]; 67 /* add non-zero cols of B into the sorted linked list lnk */ 68 ierr = PetscLLAdd(bnzj,bjj,bn,nlnk,lnk,lnkbt);CHKERRQ(ierr); 69 cnzi += nlnk; 70 } 71 72 /* If free space is not available, make more free space */ 73 /* Double the amount of total space in the list */ 74 if (current_space->local_remaining<cnzi) { 75 ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 76 nspacedouble++; 77 } 78 79 /* Copy data into free space, then initialize lnk */ 80 ierr = PetscLLClean(bn,bn,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 81 current_space->array += cnzi; 82 current_space->local_used += cnzi; 83 current_space->local_remaining -= cnzi; 84 85 ci[i+1] = ci[i] + cnzi; 86 } 87 88 /* Column indices are in the list of free space */ 89 /* Allocate space for cj, initialize cj, and */ 90 /* destroy list of free space and other temporary array(s) */ 91 ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 92 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 93 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 94 95 /* Allocate space for ca */ 96 ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 97 ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr); 98 99 /* put together the new symbolic matrix */ 100 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bn,ci,cj,ca,C);CHKERRQ(ierr); 101 102 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 103 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 104 c = (Mat_SeqAIJ *)((*C)->data); 105 c->free_a = PETSC_TRUE; 106 c->free_ij = PETSC_TRUE; 107 c->nonew = 0; 108 109 /* set MatInfo */ 110 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5; 111 if (afill < 1.0) afill = 1.0; 112 c->maxnz = ci[am]; 113 c->nz = ci[am]; 114 (*C)->info.mallocs = nspacedouble; 115 (*C)->info.fill_ratio_given = fill; 116 (*C)->info.fill_ratio_needed = afill; 117 118 #if defined(PETSC_USE_INFO) 119 if (ci[am]) { 120 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 121 ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 122 } else { 123 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 124 } 125 #endif 126 PetscFunctionReturn(0); 127 } 128 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqAIJ" 132 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 133 { 134 PetscErrorCode ierr; 135 PetscLogDouble flops=0.0; 136 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 137 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 138 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; 139 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j; 140 PetscInt am=A->rmap->N,cm=C->rmap->N; 141 PetscInt i,j,k,anzi,bnzi,cnzi,brow,nextb; 142 MatScalar *aa=a->a,*ba=b->a,*baj,*ca=c->a; 143 144 PetscFunctionBegin; 145 /* clean old values in C */ 146 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 147 /* Traverse A row-wise. */ 148 /* Build the ith row in C by summing over nonzero columns in A, */ 149 /* the rows of B corresponding to nonzeros of A. */ 150 for (i=0;i<am;i++) { 151 anzi = ai[i+1] - ai[i]; 152 for (j=0;j<anzi;j++) { 153 brow = *aj++; 154 bnzi = bi[brow+1] - bi[brow]; 155 bjj = bj + bi[brow]; 156 baj = ba + bi[brow]; 157 nextb = 0; 158 for (k=0; nextb<bnzi; k++) { 159 if (cj[k] == bjj[nextb]){ /* ccol == bcol */ 160 ca[k] += (*aa)*baj[nextb++]; 161 } 162 } 163 flops += 2*bnzi; 164 aa++; 165 } 166 cnzi = ci[i+1] - ci[i]; 167 ca += cnzi; 168 cj += cnzi; 169 } 170 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 171 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 172 173 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 174 PetscFunctionReturn(0); 175 } 176 177 #undef __FUNCT__ 178 #define __FUNCT__ "MatMatMultTranspose_SeqAIJ_SeqAIJ" 179 PetscErrorCode MatMatMultTranspose_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 180 { 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 if (scall == MAT_INITIAL_MATRIX){ 185 ierr = MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 186 } 187 ierr = MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 #undef __FUNCT__ 192 #define __FUNCT__ "PetscContainerDestroy_Mat_MatMatMultTrans" 193 PetscErrorCode PetscContainerDestroy_Mat_MatMatMultTrans(void *ptr) 194 { 195 PetscErrorCode ierr; 196 Mat_MatMatMultTrans *multtrans=(Mat_MatMatMultTrans*)ptr; 197 198 PetscFunctionBegin; 199 ierr = MatMultTransposeColoringDestroy(&multtrans->matcoloring);CHKERRQ(ierr); 200 ierr = MatDestroy(&multtrans->Bt_den);CHKERRQ(ierr); 201 ierr = MatDestroy(&multtrans->ABt_den);CHKERRQ(ierr); 202 ierr = PetscFree(multtrans);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 #undef __FUNCT__ 207 #define __FUNCT__ "MatDestroy_SeqAIJ_MatMatMultTrans" 208 PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(Mat A) 209 { 210 PetscErrorCode ierr; 211 PetscContainer container; 212 Mat_MatMatMultTrans *multtrans=PETSC_NULL; 213 214 PetscFunctionBegin; 215 ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatMultTrans",(PetscObject *)&container);CHKERRQ(ierr); 216 if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit"); 217 ierr = PetscContainerGetPointer(container,(void **)&multtrans);CHKERRQ(ierr); 218 A->ops->destroy = multtrans->destroy; 219 if (A->ops->destroy) { 220 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 221 } 222 ierr = PetscObjectCompose((PetscObject)A,"Mat_MatMatMultTrans",0);CHKERRQ(ierr); 223 PetscFunctionReturn(0); 224 } 225 226 #undef __FUNCT__ 227 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ" 228 PetscErrorCode MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 229 { 230 PetscErrorCode ierr; 231 Mat Bt; 232 PetscInt *bti,*btj; 233 Mat_MatMatMultTrans *multtrans; 234 PetscContainer container; 235 236 PetscFunctionBegin; 237 /* create symbolic Bt */ 238 ierr = MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr); 239 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,B->cmap->n,B->rmap->n,bti,btj,PETSC_NULL,&Bt);CHKERRQ(ierr); 240 241 /* get symbolic C=A*Bt */ 242 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);CHKERRQ(ierr); 243 244 /* clean up */ 245 ierr = MatDestroy(&Bt);CHKERRQ(ierr); 246 ierr = MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr); 247 248 /* create a supporting struct for reuse intermidiate dense matrices with matcoloring */ 249 ierr = PetscNew(Mat_MatMatMultTrans,&multtrans);CHKERRQ(ierr); 250 251 /* attach the supporting struct to C */ 252 ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr); 253 ierr = PetscContainerSetPointer(container,multtrans);CHKERRQ(ierr); 254 ierr = PetscContainerSetUserDestroy(container,PetscContainerDestroy_Mat_MatMatMultTrans);CHKERRQ(ierr); 255 ierr = PetscObjectCompose((PetscObject)(*C),"Mat_MatMatMultTrans",(PetscObject)container);CHKERRQ(ierr); 256 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 257 258 multtrans->usecoloring = PETSC_FALSE; 259 multtrans->destroy = (*C)->ops->destroy; 260 (*C)->ops->destroy = MatDestroy_SeqAIJ_MatMatMultTrans; 261 262 ierr = PetscOptionsGetBool(PETSC_NULL,"-matmulttrans_color",&multtrans->usecoloring,PETSC_NULL);CHKERRQ(ierr); 263 if (multtrans->usecoloring){ 264 /* Create MatMultTransposeColoring from symbolic C=A*B^T */ 265 MatMultTransposeColoring matcoloring; 266 ISColoring iscoloring; 267 Mat Bt_dense,C_dense; 268 printf("Create MatMultTransposeColoring ...\n"); 269 ierr = MatGetColoring(*C,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 270 ierr = MatMultTransposeColoringCreate(*C,iscoloring,&matcoloring);CHKERRQ(ierr); 271 multtrans->matcoloring = matcoloring; 272 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 273 printf("Create MatMultTransposeColoring is done \n"); 274 275 /* Create Bt_dense and C_dense = A*Bt_dense */ 276 ierr = MatCreate(PETSC_COMM_SELF,&Bt_dense);CHKERRQ(ierr); 277 ierr = MatSetSizes(Bt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);CHKERRQ(ierr); 278 ierr = MatSetType(Bt_dense,MATSEQDENSE);CHKERRQ(ierr); 279 ierr = MatSeqDenseSetPreallocation(Bt_dense,PETSC_NULL);CHKERRQ(ierr); 280 Bt_dense->assembled = PETSC_TRUE; 281 multtrans->Bt_den = Bt_dense; 282 283 ierr = MatCreate(PETSC_COMM_SELF,&C_dense);CHKERRQ(ierr); 284 ierr = MatSetSizes(C_dense,A->rmap->n,matcoloring->ncolors,A->rmap->n,matcoloring->ncolors);CHKERRQ(ierr); 285 ierr = MatSetType(C_dense,MATSEQDENSE);CHKERRQ(ierr); 286 ierr = MatSeqDenseSetPreallocation(C_dense,PETSC_NULL);CHKERRQ(ierr); 287 Bt_dense->assembled = PETSC_TRUE; 288 multtrans->ABt_den = C_dense; 289 } 290 291 #if defined(INEFFICIENT_ALGORITHM) 292 /* The algorithm below computes am*bm sparse inner-product - inefficient! It will be deleted later. */ 293 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 294 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 295 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*ci,*cj,*acol,*bcol; 296 PetscInt am=A->rmap->N,bm=B->rmap->N; 297 PetscInt i,j,anzi,bnzj,cnzi,nlnk,*lnk,nspacedouble=0,ka,kb,index[1]; 298 MatScalar *ca; 299 PetscBT lnkbt; 300 PetscReal afill; 301 302 /* Allocate row pointer array ci */ 303 ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 304 ci[0] = 0; 305 306 /* Create and initialize a linked list for C columns */ 307 nlnk = bm+1; 308 ierr = PetscLLCreate(bm,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr); 309 310 /* Initial FreeSpace with size fill*(nnz(A)+nnz(B)) */ 311 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr); 312 current_space = free_space; 313 314 /* Determine symbolic info for each row of the product A*B^T: */ 315 for (i=0; i<am; i++) { 316 anzi = ai[i+1] - ai[i]; 317 cnzi = 0; 318 acol = aj + ai[i]; 319 for (j=0; j<bm; j++){ 320 bnzj = bi[j+1] - bi[j]; 321 bcol= bj + bi[j]; 322 /* sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */ 323 ka = 0; kb = 0; 324 while (ka < anzi && kb < bnzj){ 325 while (acol[ka] < bcol[kb] && ka < anzi) ka++; 326 if (ka == anzi) break; 327 while (acol[ka] > bcol[kb] && kb < bnzj) kb++; 328 if (kb == bnzj) break; 329 if (acol[ka] == bcol[kb]){ /* add nonzero c(i,j) to lnk */ 330 index[0] = j; 331 ierr = PetscLLAdd(1,index,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr); 332 cnzi++; 333 break; 334 } 335 } 336 } 337 338 /* If free space is not available, make more free space */ 339 /* Double the amount of total space in the list */ 340 if (current_space->local_remaining<cnzi) { 341 ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 342 nspacedouble++; 343 } 344 345 /* Copy data into free space, then initialize lnk */ 346 ierr = PetscLLClean(bm,bm,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 347 current_space->array += cnzi; 348 current_space->local_used += cnzi; 349 current_space->local_remaining -= cnzi; 350 351 ci[i+1] = ci[i] + cnzi; 352 } 353 354 355 /* Column indices are in the list of free space. 356 Allocate array cj, copy column indices to cj, and destroy list of free space */ 357 ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 358 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 359 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 360 361 /* Allocate space for ca */ 362 ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 363 ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr); 364 365 /* put together the new symbolic matrix */ 366 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bm,ci,cj,ca,C);CHKERRQ(ierr); 367 368 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 369 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 370 c = (Mat_SeqAIJ *)((*C)->data); 371 c->free_a = PETSC_TRUE; 372 c->free_ij = PETSC_TRUE; 373 c->nonew = 0; 374 375 /* set MatInfo */ 376 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5; 377 if (afill < 1.0) afill = 1.0; 378 c->maxnz = ci[am]; 379 c->nz = ci[am]; 380 (*C)->info.mallocs = nspacedouble; 381 (*C)->info.fill_ratio_given = fill; 382 (*C)->info.fill_ratio_needed = afill; 383 384 #if defined(PETSC_USE_INFO) 385 if (ci[am]) { 386 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 387 ierr = PetscInfo1((*C),"Use MatMatMultTranspose(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 388 } else { 389 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 390 } 391 #endif 392 #endif 393 PetscFunctionReturn(0); 394 } 395 396 /* #define USE_ARRAY - for sparse dot product. Slower than !USE_ARRAY */ 397 #undef __FUNCT__ 398 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ" 399 PetscErrorCode MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 400 { 401 PetscErrorCode ierr; 402 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data; 403 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,anzi,bnzj,nexta,nextb,*acol,*bcol,brow; 404 PetscInt cm=C->rmap->n,*ci=c->i,*cj=c->j,i,j,cnzi,*ccol; 405 PetscLogDouble flops=0.0; 406 MatScalar *aa=a->a,*aval,*ba=b->a,*bval,*ca=c->a,*cval; 407 Mat_MatMatMultTrans *multtrans; 408 PetscContainer container; 409 #if defined(USE_ARRAY) 410 MatScalar *spdot; 411 #endif 412 413 PetscFunctionBegin; 414 ierr = PetscObjectQuery((PetscObject)C,"Mat_MatMatMultTrans",(PetscObject *)&container);CHKERRQ(ierr); 415 if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit"); 416 ierr = PetscContainerGetPointer(container,(void **)&multtrans);CHKERRQ(ierr); 417 if (multtrans->usecoloring){ 418 MatMultTransposeColoring matcoloring = multtrans->matcoloring; 419 Mat Bt_dense; 420 PetscInt m,n; 421 PetscLogDouble t0,tf,etime1=0.0,etime2=0.0; 422 Mat C_dense = multtrans->ABt_den; 423 424 Bt_dense = multtrans->Bt_den; 425 ierr = MatGetLocalSize(Bt_dense,&m,&n);CHKERRQ(ierr); 426 printf("Bt_dense: %d,%d\n",m,n); 427 428 /* Get Bt_dense by Apply MatMultTransposeColoring to B */ 429 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 430 ierr = MatMultTransposeColoringApply(B,Bt_dense,matcoloring);CHKERRQ(ierr); 431 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 432 etime1 += tf - t0; 433 434 /* C_dense = A*Bt_dense */ 435 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 436 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,Bt_dense,C_dense);CHKERRQ(ierr); 437 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 438 etime2 += tf - t0; 439 440 /* Recover C from C_dense */ 441 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 442 ierr = MatMultTransColoringApplyDenToSp(matcoloring,C_dense,C);CHKERRQ(ierr); 443 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 444 etime1 += tf - t0; 445 printf("etime ColoringApply: %g; MatMatMultNumeric_sp_dense: %g\n",etime1,etime2); 446 PetscFunctionReturn(0); 447 } 448 449 #if defined(USE_ARRAY) 450 /* allocate an array for implementing sparse inner-product */ 451 ierr = PetscMalloc((A->cmap->n+1)*sizeof(MatScalar),&spdot);CHKERRQ(ierr); 452 ierr = PetscMemzero(spdot,(A->cmap->n+1)*sizeof(MatScalar));CHKERRQ(ierr); 453 #endif 454 455 /* clear old values in C */ 456 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 457 458 for (i=0; i<cm; i++) { 459 anzi = ai[i+1] - ai[i]; 460 acol = aj + ai[i]; 461 aval = aa + ai[i]; 462 cnzi = ci[i+1] - ci[i]; 463 ccol = cj + ci[i]; 464 cval = ca + ci[i]; 465 for (j=0; j<cnzi; j++){ 466 brow = ccol[j]; 467 bnzj = bi[brow+1] - bi[brow]; 468 bcol = bj + bi[brow]; 469 bval = ba + bi[brow]; 470 471 /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */ 472 #if defined(USE_ARRAY) 473 /* put ba to spdot array */ 474 for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = bval[nextb]; 475 /* c(i,j)=A[i,:]*B[j,:]^T */ 476 for (nexta=0; nexta<anzi; nexta++){ 477 cval[j] += spdot[acol[nexta]]*aval[nexta]; 478 } 479 /* zero spdot array */ 480 for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = 0.0; 481 #else 482 nexta = 0; nextb = 0; 483 while (nexta<anzi && nextb<bnzj){ 484 while (acol[nexta] < bcol[nextb] && nexta < anzi) nexta++; 485 if (nexta == anzi) break; 486 while (acol[nexta] > bcol[nextb] && nextb < bnzj) nextb++; 487 if (nextb == bnzj) break; 488 if (acol[nexta] == bcol[nextb]){ 489 cval[j] += aval[nexta]*bval[nextb]; 490 nexta++; nextb++; 491 flops += 2; 492 } 493 } 494 #endif 495 } 496 } 497 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 498 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 499 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 500 #if defined(USE_ARRAY) 501 ierr = PetscFree(spdot); 502 #endif 503 PetscFunctionReturn(0); 504 } 505 506 #undef __FUNCT__ 507 #define __FUNCT__ "MatMatTransposeMult_SeqAIJ_SeqAIJ" 508 PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) { 509 PetscErrorCode ierr; 510 511 PetscFunctionBegin; 512 if (scall == MAT_INITIAL_MATRIX){ 513 ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 514 } 515 ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 516 PetscFunctionReturn(0); 517 } 518 519 #undef __FUNCT__ 520 #define __FUNCT__ "MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ" 521 PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 522 { 523 PetscErrorCode ierr; 524 Mat At; 525 PetscInt *ati,*atj; 526 527 PetscFunctionBegin; 528 /* create symbolic At */ 529 ierr = MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 530 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->cmap->n,A->rmap->n,ati,atj,PETSC_NULL,&At);CHKERRQ(ierr); 531 532 /* get symbolic C=At*B */ 533 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);CHKERRQ(ierr); 534 535 /* clean up */ 536 ierr = MatDestroy(&At);CHKERRQ(ierr); 537 ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 538 PetscFunctionReturn(0); 539 } 540 541 #undef __FUNCT__ 542 #define __FUNCT__ "MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ" 543 PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 544 { 545 PetscErrorCode ierr; 546 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data; 547 PetscInt am=A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb; 548 PetscInt cm=C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k; 549 PetscLogDouble flops=0.0; 550 MatScalar *aa=a->a,*ba,*ca=c->a,*caj; 551 552 PetscFunctionBegin; 553 /* clear old values in C */ 554 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 555 556 /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */ 557 for (i=0;i<am;i++) { 558 bj = b->j + bi[i]; 559 ba = b->a + bi[i]; 560 bnzi = bi[i+1] - bi[i]; 561 anzi = ai[i+1] - ai[i]; 562 for (j=0; j<anzi; j++) { 563 nextb = 0; 564 crow = *aj++; 565 cjj = cj + ci[crow]; 566 caj = ca + ci[crow]; 567 /* perform sparse axpy operation. Note cjj includes bj. */ 568 for (k=0; nextb<bnzi; k++) { 569 if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */ 570 caj[k] += (*aa)*(*(ba+nextb)); 571 nextb++; 572 } 573 } 574 flops += 2*bnzi; 575 aa++; 576 } 577 } 578 579 /* Assemble the final matrix and clean up */ 580 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 581 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 582 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 583 PetscFunctionReturn(0); 584 } 585 586 EXTERN_C_BEGIN 587 #undef __FUNCT__ 588 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqDense" 589 PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 590 { 591 PetscErrorCode ierr; 592 593 PetscFunctionBegin; 594 if (scall == MAT_INITIAL_MATRIX){ 595 ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);CHKERRQ(ierr); 596 } 597 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);CHKERRQ(ierr); 598 PetscFunctionReturn(0); 599 } 600 EXTERN_C_END 601 602 #undef __FUNCT__ 603 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqDense" 604 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 605 { 606 PetscErrorCode ierr; 607 608 PetscFunctionBegin; 609 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);CHKERRQ(ierr); 610 PetscFunctionReturn(0); 611 } 612 613 #undef __FUNCT__ 614 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqDense" 615 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 616 { 617 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 618 PetscErrorCode ierr; 619 PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 620 MatScalar *aa; 621 PetscInt cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n; 622 PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam; 623 624 PetscFunctionBegin; 625 if (!cm || !cn) PetscFunctionReturn(0); 626 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); 627 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); 628 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); 629 ierr = MatGetArray(B,&b);CHKERRQ(ierr); 630 ierr = MatGetArray(C,&c);CHKERRQ(ierr); 631 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 632 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 633 colam = col*am; 634 for (i=0; i<am; i++) { /* over rows of C in those columns */ 635 r1 = r2 = r3 = r4 = 0.0; 636 n = a->i[i+1] - a->i[i]; 637 aj = a->j + a->i[i]; 638 aa = a->a + a->i[i]; 639 for (j=0; j<n; j++) { 640 r1 += (*aa)*b1[*aj]; 641 r2 += (*aa)*b2[*aj]; 642 r3 += (*aa)*b3[*aj]; 643 r4 += (*aa++)*b4[*aj++]; 644 } 645 c[colam + i] = r1; 646 c[colam + am + i] = r2; 647 c[colam + am2 + i] = r3; 648 c[colam + am3 + i] = r4; 649 } 650 b1 += bm4; 651 b2 += bm4; 652 b3 += bm4; 653 b4 += bm4; 654 } 655 for (;col<cn; col++){ /* over extra columns of C */ 656 for (i=0; i<am; i++) { /* over rows of C in those columns */ 657 r1 = 0.0; 658 n = a->i[i+1] - a->i[i]; 659 aj = a->j + a->i[i]; 660 aa = a->a + a->i[i]; 661 662 for (j=0; j<n; j++) { 663 r1 += (*aa++)*b1[*aj++]; 664 } 665 c[col*am + i] = r1; 666 } 667 b1 += bm; 668 } 669 ierr = PetscLogFlops(cn*(2.0*a->nz));CHKERRQ(ierr); 670 ierr = MatRestoreArray(B,&b);CHKERRQ(ierr); 671 ierr = MatRestoreArray(C,&c);CHKERRQ(ierr); 672 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 673 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 674 PetscFunctionReturn(0); 675 } 676 677 /* 678 Note very similar to MatMult_SeqAIJ(), should generate both codes from same base 679 */ 680 #undef __FUNCT__ 681 #define __FUNCT__ "MatMatMultNumericAdd_SeqAIJ_SeqDense" 682 PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 683 { 684 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 685 PetscErrorCode ierr; 686 PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 687 MatScalar *aa; 688 PetscInt cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n,*ii,arm; 689 PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam,*ridx; 690 691 PetscFunctionBegin; 692 if (!cm || !cn) PetscFunctionReturn(0); 693 ierr = MatGetArray(B,&b);CHKERRQ(ierr); 694 ierr = MatGetArray(C,&c);CHKERRQ(ierr); 695 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 696 697 if (a->compressedrow.use){ /* use compressed row format */ 698 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 699 colam = col*am; 700 arm = a->compressedrow.nrows; 701 ii = a->compressedrow.i; 702 ridx = a->compressedrow.rindex; 703 for (i=0; i<arm; i++) { /* over rows of C in those columns */ 704 r1 = r2 = r3 = r4 = 0.0; 705 n = ii[i+1] - ii[i]; 706 aj = a->j + ii[i]; 707 aa = a->a + ii[i]; 708 for (j=0; j<n; j++) { 709 r1 += (*aa)*b1[*aj]; 710 r2 += (*aa)*b2[*aj]; 711 r3 += (*aa)*b3[*aj]; 712 r4 += (*aa++)*b4[*aj++]; 713 } 714 c[colam + ridx[i]] += r1; 715 c[colam + am + ridx[i]] += r2; 716 c[colam + am2 + ridx[i]] += r3; 717 c[colam + am3 + ridx[i]] += r4; 718 } 719 b1 += bm4; 720 b2 += bm4; 721 b3 += bm4; 722 b4 += bm4; 723 } 724 for (;col<cn; col++){ /* over extra columns of C */ 725 colam = col*am; 726 arm = a->compressedrow.nrows; 727 ii = a->compressedrow.i; 728 ridx = a->compressedrow.rindex; 729 for (i=0; i<arm; i++) { /* over rows of C in those columns */ 730 r1 = 0.0; 731 n = ii[i+1] - ii[i]; 732 aj = a->j + ii[i]; 733 aa = a->a + ii[i]; 734 735 for (j=0; j<n; j++) { 736 r1 += (*aa++)*b1[*aj++]; 737 } 738 c[col*am + ridx[i]] += r1; 739 } 740 b1 += bm; 741 } 742 } else { 743 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 744 colam = col*am; 745 for (i=0; i<am; i++) { /* over rows of C in those columns */ 746 r1 = r2 = r3 = r4 = 0.0; 747 n = a->i[i+1] - a->i[i]; 748 aj = a->j + a->i[i]; 749 aa = a->a + a->i[i]; 750 for (j=0; j<n; j++) { 751 r1 += (*aa)*b1[*aj]; 752 r2 += (*aa)*b2[*aj]; 753 r3 += (*aa)*b3[*aj]; 754 r4 += (*aa++)*b4[*aj++]; 755 } 756 c[colam + i] += r1; 757 c[colam + am + i] += r2; 758 c[colam + am2 + i] += r3; 759 c[colam + am3 + i] += r4; 760 } 761 b1 += bm4; 762 b2 += bm4; 763 b3 += bm4; 764 b4 += bm4; 765 } 766 for (;col<cn; col++){ /* over extra columns of C */ 767 for (i=0; i<am; i++) { /* over rows of C in those columns */ 768 r1 = 0.0; 769 n = a->i[i+1] - a->i[i]; 770 aj = a->j + a->i[i]; 771 aa = a->a + a->i[i]; 772 773 for (j=0; j<n; j++) { 774 r1 += (*aa++)*b1[*aj++]; 775 } 776 c[col*am + i] += r1; 777 } 778 b1 += bm; 779 } 780 } 781 ierr = PetscLogFlops(cn*2.0*a->nz);CHKERRQ(ierr); 782 ierr = MatRestoreArray(B,&b);CHKERRQ(ierr); 783 ierr = MatRestoreArray(C,&c);CHKERRQ(ierr); 784 PetscFunctionReturn(0); 785 } 786 787 #undef __FUNCT__ 788 #define __FUNCT__ "MatMultTransposeColoringApply_SeqAIJ" 789 PetscErrorCode MatMultTransposeColoringApply_SeqAIJ(Mat B,Mat Btdense,MatMultTransposeColoring coloring) 790 { 791 PetscErrorCode ierr; 792 Mat_SeqAIJ *a = (Mat_SeqAIJ*)B->data; 793 Mat_SeqDense *atdense = (Mat_SeqDense*)Btdense->data; 794 PetscInt m=Btdense->rmap->n,n=Btdense->cmap->n,j,k,l,col,anz,*atcol,brow,bcol; 795 MatScalar *atval,*bval; 796 797 PetscFunctionBegin; 798 ierr = PetscMemzero(atdense->v,(m*n)*sizeof(MatScalar));CHKERRQ(ierr); 799 for (k=0; k<coloring->ncolors; k++) { 800 for (l=0; l<coloring->ncolumns[k]; l++) { /* insert a row of B to a column of Btdense */ 801 col = coloring->columns[k][l]; /* =row of B */ 802 anz = a->i[col+1] - a->i[col]; 803 for (j=0; j<anz; j++){ 804 atcol = a->j + a->i[col]; 805 atval = a->a + a->i[col]; 806 bval = atdense->v; 807 brow = atcol[j]; 808 bcol = k; 809 bval[bcol*m+brow] = atval[j]; 810 } 811 } 812 } 813 PetscFunctionReturn(0); 814 } 815 816 #undef __FUNCT__ 817 #define __FUNCT__ "MatMultTransColoringApplyDenToSp_SeqAIJ" 818 PetscErrorCode MatMultTransColoringApplyDenToSp_SeqAIJ(MatMultTransposeColoring matcoloring,Mat Cden,Mat Csp) 819 { 820 PetscErrorCode ierr; 821 PetscInt k,l,row,col,m; 822 PetscScalar *ca,*cval; 823 824 PetscFunctionBegin; 825 ierr = MatGetLocalSize(Csp,&m,PETSC_NULL);CHKERRQ(ierr); 826 ierr = MatGetArray(Cden,&ca);CHKERRQ(ierr); 827 cval = ca; 828 for (k=0; k<matcoloring->ncolors; k++) { 829 for (l=0; l<matcoloring->nrows[k]; l++){ 830 row = matcoloring->rows[k][l]; /* local row index */ 831 col = matcoloring->columnsforrow[k][l]; /* global column index */ 832 ierr = MatSetValues(Csp,1,&row,1,&col,cval+row,INSERT_VALUES);CHKERRQ(ierr); 833 } 834 cval += m; 835 } 836 ierr = MatRestoreArray(Cden,&ca);CHKERRQ(ierr); 837 ierr = MatAssemblyBegin(Csp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 838 ierr = MatAssemblyEnd(Csp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 839 PetscFunctionReturn(0); 840 } 841 842 #undef __FUNCT__ 843 #define __FUNCT__ "MatMultTransposeColoringCreate_SeqAIJ" 844 PetscErrorCode MatMultTransposeColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatMultTransposeColoring c) 845 { 846 PetscErrorCode ierr; 847 PetscInt i,n,nrows,N,j,k,m,*rows,*ci,*cj,ncols,col; 848 const PetscInt *is; 849 PetscInt nis = iscoloring->n,*rowhit,*columnsforrow,bs = 1; 850 IS *isa; 851 PetscBool done; 852 PetscBool flg1,flg2; 853 854 PetscFunctionBegin; 855 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();"); 856 857 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 858 /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */ 859 ierr = PetscTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);CHKERRQ(ierr); 860 ierr = PetscTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);CHKERRQ(ierr); 861 if (flg1 || flg2) { 862 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 863 } 864 865 N = mat->cmap->N/bs; 866 c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */ 867 c->N = mat->cmap->N/bs; 868 c->m = mat->rmap->N/bs; 869 c->rstart = 0; 870 871 c->ncolors = nis; 872 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 873 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 874 ierr = PetscMalloc(c->m*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 875 ierr = PetscMalloc(c->m*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 876 ierr = PetscMalloc(c->m*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 877 878 ierr = MatGetColumnIJ(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 879 if (!done) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"MatGetColumnIJ() not supported for matrix type %s",((PetscObject)mat)->type_name); 880 881 ierr = PetscMalloc((c->m+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 882 ierr = PetscMalloc((c->m+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr); 883 884 for (i=0; i<nis; i++) { 885 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 886 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 887 c->ncolumns[i] = n; 888 if (n) { 889 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 890 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 891 } else { 892 c->columns[i] = 0; 893 } 894 895 /* fast, crude version requires O(N*N) work */ 896 ierr = PetscMemzero(rowhit,c->m*sizeof(PetscInt));CHKERRQ(ierr); 897 /* loop over columns*/ 898 for (j=0; j<n; j++) { 899 col = is[j]; 900 rows = cj + ci[col]; 901 m = ci[col+1] - ci[col]; 902 /* loop over columns marking them in rowhit */ 903 for (k=0; k<m; k++) { 904 rowhit[*rows++] = col + 1; 905 } 906 } 907 /* count the number of hits */ 908 nrows = 0; 909 for (j=0; j<c->m; j++) { 910 if (rowhit[j]) nrows++; 911 } 912 c->nrows[i] = nrows; 913 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 914 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 915 nrows = 0; 916 for (j=0; j<c->m; j++) { 917 if (rowhit[j]) { 918 c->rows[i][nrows] = j; 919 c->columnsforrow[i][nrows] = rowhit[j] - 1; 920 nrows++; 921 } 922 } 923 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 924 } 925 ierr = MatRestoreColumnIJ(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 926 927 ierr = PetscFree(rowhit);CHKERRQ(ierr); 928 ierr = PetscFree(columnsforrow);CHKERRQ(ierr); 929 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 930 PetscFunctionReturn(0); 931 } 932