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 /* This routine is not used. Should be removed! */ 178 #undef __FUNCT__ 179 #define __FUNCT__ "MatMatTransposeMult_SeqAIJ_SeqAIJ" 180 PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 181 { 182 PetscErrorCode ierr; 183 184 PetscFunctionBegin; 185 if (scall == MAT_INITIAL_MATRIX){ 186 ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 187 } 188 ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191 192 #undef __FUNCT__ 193 #define __FUNCT__ "PetscContainerDestroy_Mat_MatMatMultTrans" 194 PetscErrorCode PetscContainerDestroy_Mat_MatMatMultTrans(void *ptr) 195 { 196 PetscErrorCode ierr; 197 Mat_MatMatMultTrans *multtrans=(Mat_MatMatMultTrans*)ptr; 198 199 PetscFunctionBegin; 200 ierr = MatTransposeColoringDestroy(&multtrans->matcoloring);CHKERRQ(ierr); 201 ierr = MatDestroy(&multtrans->Bt_den);CHKERRQ(ierr); 202 ierr = MatDestroy(&multtrans->ABt_den);CHKERRQ(ierr); 203 ierr = PetscFree(multtrans);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "MatDestroy_SeqAIJ_MatMatMultTrans" 209 PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(Mat A) 210 { 211 PetscErrorCode ierr; 212 PetscContainer container; 213 Mat_MatMatMultTrans *multtrans=PETSC_NULL; 214 215 PetscFunctionBegin; 216 ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatMultTrans",(PetscObject *)&container);CHKERRQ(ierr); 217 if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit"); 218 ierr = PetscContainerGetPointer(container,(void **)&multtrans);CHKERRQ(ierr); 219 A->ops->destroy = multtrans->destroy; 220 if (A->ops->destroy) { 221 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 222 } 223 ierr = PetscObjectCompose((PetscObject)A,"Mat_MatMatMultTrans",0);CHKERRQ(ierr); 224 PetscFunctionReturn(0); 225 } 226 227 #undef __FUNCT__ 228 #define __FUNCT__ "MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ" 229 PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 230 { 231 PetscErrorCode ierr; 232 Mat Bt; 233 PetscInt *bti,*btj; 234 Mat_MatMatMultTrans *multtrans; 235 PetscContainer container; 236 237 PetscLogDouble t0,tf,etime2=0.0; 238 239 PetscFunctionBegin; 240 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 241 /* create symbolic Bt */ 242 ierr = MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr); 243 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,B->cmap->n,B->rmap->n,bti,btj,PETSC_NULL,&Bt);CHKERRQ(ierr); 244 245 /* get symbolic C=A*Bt */ 246 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);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 = PetscGetTime(&tf);CHKERRQ(ierr); 263 etime2 += tf - t0; 264 265 ierr = PetscOptionsGetBool(PETSC_NULL,"-matmulttrans_color",&multtrans->usecoloring,PETSC_NULL);CHKERRQ(ierr); 266 if (multtrans->usecoloring){ 267 /* Create MatTransposeColoring from symbolic C=A*B^T */ 268 MatTransposeColoring matcoloring; 269 ISColoring iscoloring; 270 Mat Bt_dense,C_dense; 271 PetscLogDouble etime0=0.0,etime01=0.0,etime1=0.0; 272 273 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 274 ierr = MatGetColoring(*C,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 275 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 276 etime0 += tf - t0; 277 278 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 279 ierr = MatTransposeColoringCreate(*C,iscoloring,&matcoloring);CHKERRQ(ierr); 280 multtrans->matcoloring = matcoloring; 281 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 282 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 283 etime01 += tf - t0; 284 285 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 286 /* Create Bt_dense and C_dense = A*Bt_dense */ 287 ierr = MatCreate(PETSC_COMM_SELF,&Bt_dense);CHKERRQ(ierr); 288 ierr = MatSetSizes(Bt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);CHKERRQ(ierr); 289 ierr = MatSetType(Bt_dense,MATSEQDENSE);CHKERRQ(ierr); 290 ierr = MatSeqDenseSetPreallocation(Bt_dense,PETSC_NULL);CHKERRQ(ierr); 291 Bt_dense->assembled = PETSC_TRUE; 292 multtrans->Bt_den = Bt_dense; 293 294 ierr = MatCreate(PETSC_COMM_SELF,&C_dense);CHKERRQ(ierr); 295 ierr = MatSetSizes(C_dense,A->rmap->n,matcoloring->ncolors,A->rmap->n,matcoloring->ncolors);CHKERRQ(ierr); 296 ierr = MatSetType(C_dense,MATSEQDENSE);CHKERRQ(ierr); 297 ierr = MatSeqDenseSetPreallocation(C_dense,PETSC_NULL);CHKERRQ(ierr); 298 Bt_dense->assembled = PETSC_TRUE; 299 multtrans->ABt_den = C_dense; 300 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 301 etime1 += tf - t0; 302 printf("Sym = GetColor %g + ColorCreate %g + MatDenseCreate %g + non-colorSym %g = %g\n",etime0,etime01,etime1,etime2,etime0+etime01+etime1+etime2); 303 } 304 /* clean up */ 305 ierr = MatDestroy(&Bt);CHKERRQ(ierr); 306 ierr = MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr); 307 308 #if defined(INEFFICIENT_ALGORITHM) 309 /* The algorithm below computes am*bm sparse inner-product - inefficient! It will be deleted later. */ 310 PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL; 311 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c; 312 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*ci,*cj,*acol,*bcol; 313 PetscInt am=A->rmap->N,bm=B->rmap->N; 314 PetscInt i,j,anzi,bnzj,cnzi,nlnk,*lnk,nspacedouble=0,ka,kb,index[1]; 315 MatScalar *ca; 316 PetscBT lnkbt; 317 PetscReal afill; 318 319 /* Allocate row pointer array ci */ 320 ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr); 321 ci[0] = 0; 322 323 /* Create and initialize a linked list for C columns */ 324 nlnk = bm+1; 325 ierr = PetscLLCreate(bm,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr); 326 327 /* Initial FreeSpace with size fill*(nnz(A)+nnz(B)) */ 328 ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr); 329 current_space = free_space; 330 331 /* Determine symbolic info for each row of the product A*B^T: */ 332 for (i=0; i<am; i++) { 333 anzi = ai[i+1] - ai[i]; 334 cnzi = 0; 335 acol = aj + ai[i]; 336 for (j=0; j<bm; j++){ 337 bnzj = bi[j+1] - bi[j]; 338 bcol= bj + bi[j]; 339 /* sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */ 340 ka = 0; kb = 0; 341 while (ka < anzi && kb < bnzj){ 342 while (acol[ka] < bcol[kb] && ka < anzi) ka++; 343 if (ka == anzi) break; 344 while (acol[ka] > bcol[kb] && kb < bnzj) kb++; 345 if (kb == bnzj) break; 346 if (acol[ka] == bcol[kb]){ /* add nonzero c(i,j) to lnk */ 347 index[0] = j; 348 ierr = PetscLLAdd(1,index,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr); 349 cnzi++; 350 break; 351 } 352 } 353 } 354 355 /* If free space is not available, make more free space */ 356 /* Double the amount of total space in the list */ 357 if (current_space->local_remaining<cnzi) { 358 ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,¤t_space);CHKERRQ(ierr); 359 nspacedouble++; 360 } 361 362 /* Copy data into free space, then initialize lnk */ 363 ierr = PetscLLClean(bm,bm,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 364 current_space->array += cnzi; 365 current_space->local_used += cnzi; 366 current_space->local_remaining -= cnzi; 367 368 ci[i+1] = ci[i] + cnzi; 369 } 370 371 372 /* Column indices are in the list of free space. 373 Allocate array cj, copy column indices to cj, and destroy list of free space */ 374 ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr); 375 ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr); 376 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 377 378 /* Allocate space for ca */ 379 ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr); 380 ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr); 381 382 /* put together the new symbolic matrix */ 383 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bm,ci,cj,ca,C);CHKERRQ(ierr); 384 385 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 386 /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ 387 c = (Mat_SeqAIJ *)((*C)->data); 388 c->free_a = PETSC_TRUE; 389 c->free_ij = PETSC_TRUE; 390 c->nonew = 0; 391 392 /* set MatInfo */ 393 afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5; 394 if (afill < 1.0) afill = 1.0; 395 c->maxnz = ci[am]; 396 c->nz = ci[am]; 397 (*C)->info.mallocs = nspacedouble; 398 (*C)->info.fill_ratio_given = fill; 399 (*C)->info.fill_ratio_needed = afill; 400 401 #if defined(PETSC_USE_INFO) 402 if (ci[am]) { 403 ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr); 404 ierr = PetscInfo1((*C),"Use MatMatTransposeMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr); 405 } else { 406 ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr); 407 } 408 #endif 409 #endif 410 PetscFunctionReturn(0); 411 } 412 413 /* #define USE_ARRAY - for sparse dot product. Slower than !USE_ARRAY */ 414 #undef __FUNCT__ 415 #define __FUNCT__ "MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ" 416 PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 417 { 418 PetscErrorCode ierr; 419 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data; 420 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,anzi,bnzj,nexta,nextb,*acol,*bcol,brow; 421 PetscInt cm=C->rmap->n,*ci=c->i,*cj=c->j,i,j,cnzi,*ccol; 422 PetscLogDouble flops=0.0; 423 MatScalar *aa=a->a,*aval,*ba=b->a,*bval,*ca=c->a,*cval; 424 Mat_MatMatMultTrans *multtrans; 425 PetscContainer container; 426 #if defined(USE_ARRAY) 427 MatScalar *spdot; 428 #endif 429 430 PetscFunctionBegin; 431 ierr = PetscObjectQuery((PetscObject)C,"Mat_MatMatMultTrans",(PetscObject *)&container);CHKERRQ(ierr); 432 if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit"); 433 ierr = PetscContainerGetPointer(container,(void **)&multtrans);CHKERRQ(ierr); 434 if (multtrans->usecoloring){ 435 MatTransposeColoring matcoloring = multtrans->matcoloring; 436 Mat Bt_dense; 437 PetscInt m,n; 438 PetscLogDouble t0,tf,etime0=0.0,etime1=0.0,etime2=0.0; 439 Mat C_dense = multtrans->ABt_den; 440 441 Bt_dense = multtrans->Bt_den; 442 ierr = MatGetLocalSize(Bt_dense,&m,&n);CHKERRQ(ierr); 443 printf("Bt_dense: %d,%d\n",m,n); 444 445 /* Get Bt_dense by Apply MatTransposeColoring to B */ 446 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 447 ierr = MatTransColoringApplySpToDen(matcoloring,B,Bt_dense);CHKERRQ(ierr); 448 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 449 etime0 += tf - t0; 450 451 /* C_dense = A*Bt_dense */ 452 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 453 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,Bt_dense,C_dense);CHKERRQ(ierr); 454 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 455 etime2 += tf - t0; 456 457 /* Recover C from C_dense */ 458 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 459 ierr = MatTransColoringApplyDenToSp(matcoloring,C_dense,C);CHKERRQ(ierr); 460 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 461 etime1 += tf - t0; 462 printf("Num = ColoringApply: %g %g + Mult_sp_dense %g = %g\n",etime0,etime1,etime2,etime0+etime1+etime2); 463 PetscFunctionReturn(0); 464 } 465 466 #if defined(USE_ARRAY) 467 /* allocate an array for implementing sparse inner-product */ 468 ierr = PetscMalloc((A->cmap->n+1)*sizeof(MatScalar),&spdot);CHKERRQ(ierr); 469 ierr = PetscMemzero(spdot,(A->cmap->n+1)*sizeof(MatScalar));CHKERRQ(ierr); 470 #endif 471 472 /* clear old values in C */ 473 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 474 475 for (i=0; i<cm; i++) { 476 anzi = ai[i+1] - ai[i]; 477 acol = aj + ai[i]; 478 aval = aa + ai[i]; 479 cnzi = ci[i+1] - ci[i]; 480 ccol = cj + ci[i]; 481 cval = ca + ci[i]; 482 for (j=0; j<cnzi; j++){ 483 brow = ccol[j]; 484 bnzj = bi[brow+1] - bi[brow]; 485 bcol = bj + bi[brow]; 486 bval = ba + bi[brow]; 487 488 /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */ 489 #if defined(USE_ARRAY) 490 /* put ba to spdot array */ 491 for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = bval[nextb]; 492 /* c(i,j)=A[i,:]*B[j,:]^T */ 493 for (nexta=0; nexta<anzi; nexta++){ 494 cval[j] += spdot[acol[nexta]]*aval[nexta]; 495 } 496 /* zero spdot array */ 497 for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = 0.0; 498 #else 499 nexta = 0; nextb = 0; 500 while (nexta<anzi && nextb<bnzj){ 501 while (acol[nexta] < bcol[nextb] && nexta < anzi) nexta++; 502 if (nexta == anzi) break; 503 while (acol[nexta] > bcol[nextb] && nextb < bnzj) nextb++; 504 if (nextb == bnzj) break; 505 if (acol[nexta] == bcol[nextb]){ 506 cval[j] += aval[nexta]*bval[nextb]; 507 nexta++; nextb++; 508 flops += 2; 509 } 510 } 511 #endif 512 } 513 } 514 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 515 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 516 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 517 #if defined(USE_ARRAY) 518 ierr = PetscFree(spdot); 519 #endif 520 PetscFunctionReturn(0); 521 } 522 523 #undef __FUNCT__ 524 #define __FUNCT__ "MatTransposeMatMult_SeqAIJ_SeqAIJ" 525 PetscErrorCode MatTransposeMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) { 526 PetscErrorCode ierr; 527 528 PetscFunctionBegin; 529 if (scall == MAT_INITIAL_MATRIX){ 530 ierr = MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 531 } 532 ierr = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr); 533 PetscFunctionReturn(0); 534 } 535 536 #undef __FUNCT__ 537 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ" 538 PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 539 { 540 PetscErrorCode ierr; 541 Mat At; 542 PetscInt *ati,*atj; 543 544 PetscFunctionBegin; 545 /* create symbolic At */ 546 ierr = MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 547 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->cmap->n,A->rmap->n,ati,atj,PETSC_NULL,&At);CHKERRQ(ierr); 548 549 /* get symbolic C=At*B */ 550 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);CHKERRQ(ierr); 551 552 /* clean up */ 553 ierr = MatDestroy(&At);CHKERRQ(ierr); 554 ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr); 555 PetscFunctionReturn(0); 556 } 557 558 #undef __FUNCT__ 559 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ" 560 PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C) 561 { 562 PetscErrorCode ierr; 563 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data; 564 PetscInt am=A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb; 565 PetscInt cm=C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k; 566 PetscLogDouble flops=0.0; 567 MatScalar *aa=a->a,*ba,*ca=c->a,*caj; 568 569 PetscFunctionBegin; 570 /* clear old values in C */ 571 ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr); 572 573 /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */ 574 for (i=0;i<am;i++) { 575 bj = b->j + bi[i]; 576 ba = b->a + bi[i]; 577 bnzi = bi[i+1] - bi[i]; 578 anzi = ai[i+1] - ai[i]; 579 for (j=0; j<anzi; j++) { 580 nextb = 0; 581 crow = *aj++; 582 cjj = cj + ci[crow]; 583 caj = ca + ci[crow]; 584 /* perform sparse axpy operation. Note cjj includes bj. */ 585 for (k=0; nextb<bnzi; k++) { 586 if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */ 587 caj[k] += (*aa)*(*(ba+nextb)); 588 nextb++; 589 } 590 } 591 flops += 2*bnzi; 592 aa++; 593 } 594 } 595 596 /* Assemble the final matrix and clean up */ 597 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 598 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 599 ierr = PetscLogFlops(flops);CHKERRQ(ierr); 600 PetscFunctionReturn(0); 601 } 602 603 EXTERN_C_BEGIN 604 #undef __FUNCT__ 605 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqDense" 606 PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 607 { 608 PetscErrorCode ierr; 609 610 PetscFunctionBegin; 611 if (scall == MAT_INITIAL_MATRIX){ 612 ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);CHKERRQ(ierr); 613 } 614 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);CHKERRQ(ierr); 615 PetscFunctionReturn(0); 616 } 617 EXTERN_C_END 618 619 #undef __FUNCT__ 620 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqDense" 621 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 622 { 623 PetscErrorCode ierr; 624 625 PetscFunctionBegin; 626 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);CHKERRQ(ierr); 627 PetscFunctionReturn(0); 628 } 629 630 #undef __FUNCT__ 631 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqDense" 632 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 633 { 634 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 635 PetscErrorCode ierr; 636 PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 637 MatScalar *aa; 638 PetscInt cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n; 639 PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam; 640 641 PetscFunctionBegin; 642 if (!cm || !cn) PetscFunctionReturn(0); 643 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); 644 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); 645 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); 646 ierr = MatGetArray(B,&b);CHKERRQ(ierr); 647 ierr = MatGetArray(C,&c);CHKERRQ(ierr); 648 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 649 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 650 colam = col*am; 651 for (i=0; i<am; i++) { /* over rows of C in those columns */ 652 r1 = r2 = r3 = r4 = 0.0; 653 n = a->i[i+1] - a->i[i]; 654 aj = a->j + a->i[i]; 655 aa = a->a + a->i[i]; 656 for (j=0; j<n; j++) { 657 r1 += (*aa)*b1[*aj]; 658 r2 += (*aa)*b2[*aj]; 659 r3 += (*aa)*b3[*aj]; 660 r4 += (*aa++)*b4[*aj++]; 661 } 662 c[colam + i] = r1; 663 c[colam + am + i] = r2; 664 c[colam + am2 + i] = r3; 665 c[colam + am3 + i] = r4; 666 } 667 b1 += bm4; 668 b2 += bm4; 669 b3 += bm4; 670 b4 += bm4; 671 } 672 for (;col<cn; col++){ /* over extra columns of C */ 673 for (i=0; i<am; i++) { /* over rows of C in those columns */ 674 r1 = 0.0; 675 n = a->i[i+1] - a->i[i]; 676 aj = a->j + a->i[i]; 677 aa = a->a + a->i[i]; 678 679 for (j=0; j<n; j++) { 680 r1 += (*aa++)*b1[*aj++]; 681 } 682 c[col*am + i] = r1; 683 } 684 b1 += bm; 685 } 686 ierr = PetscLogFlops(cn*(2.0*a->nz));CHKERRQ(ierr); 687 ierr = MatRestoreArray(B,&b);CHKERRQ(ierr); 688 ierr = MatRestoreArray(C,&c);CHKERRQ(ierr); 689 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 690 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 691 PetscFunctionReturn(0); 692 } 693 694 /* 695 Note very similar to MatMult_SeqAIJ(), should generate both codes from same base 696 */ 697 #undef __FUNCT__ 698 #define __FUNCT__ "MatMatMultNumericAdd_SeqAIJ_SeqDense" 699 PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C) 700 { 701 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 702 PetscErrorCode ierr; 703 PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 704 MatScalar *aa; 705 PetscInt cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n,*ii,arm; 706 PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam,*ridx; 707 708 PetscFunctionBegin; 709 if (!cm || !cn) PetscFunctionReturn(0); 710 ierr = MatGetArray(B,&b);CHKERRQ(ierr); 711 ierr = MatGetArray(C,&c);CHKERRQ(ierr); 712 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 713 714 if (a->compressedrow.use){ /* use compressed row format */ 715 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 716 colam = col*am; 717 arm = a->compressedrow.nrows; 718 ii = a->compressedrow.i; 719 ridx = a->compressedrow.rindex; 720 for (i=0; i<arm; i++) { /* over rows of C in those columns */ 721 r1 = r2 = r3 = r4 = 0.0; 722 n = ii[i+1] - ii[i]; 723 aj = a->j + ii[i]; 724 aa = a->a + ii[i]; 725 for (j=0; j<n; j++) { 726 r1 += (*aa)*b1[*aj]; 727 r2 += (*aa)*b2[*aj]; 728 r3 += (*aa)*b3[*aj]; 729 r4 += (*aa++)*b4[*aj++]; 730 } 731 c[colam + ridx[i]] += r1; 732 c[colam + am + ridx[i]] += r2; 733 c[colam + am2 + ridx[i]] += r3; 734 c[colam + am3 + ridx[i]] += r4; 735 } 736 b1 += bm4; 737 b2 += bm4; 738 b3 += bm4; 739 b4 += bm4; 740 } 741 for (;col<cn; col++){ /* over extra columns of C */ 742 colam = col*am; 743 arm = a->compressedrow.nrows; 744 ii = a->compressedrow.i; 745 ridx = a->compressedrow.rindex; 746 for (i=0; i<arm; i++) { /* over rows of C in those columns */ 747 r1 = 0.0; 748 n = ii[i+1] - ii[i]; 749 aj = a->j + ii[i]; 750 aa = a->a + ii[i]; 751 752 for (j=0; j<n; j++) { 753 r1 += (*aa++)*b1[*aj++]; 754 } 755 c[col*am + ridx[i]] += r1; 756 } 757 b1 += bm; 758 } 759 } else { 760 for (col=0; col<cn-4; col += 4){ /* over columns of C */ 761 colam = col*am; 762 for (i=0; i<am; i++) { /* over rows of C in those columns */ 763 r1 = r2 = r3 = r4 = 0.0; 764 n = a->i[i+1] - a->i[i]; 765 aj = a->j + a->i[i]; 766 aa = a->a + a->i[i]; 767 for (j=0; j<n; j++) { 768 r1 += (*aa)*b1[*aj]; 769 r2 += (*aa)*b2[*aj]; 770 r3 += (*aa)*b3[*aj]; 771 r4 += (*aa++)*b4[*aj++]; 772 } 773 c[colam + i] += r1; 774 c[colam + am + i] += r2; 775 c[colam + am2 + i] += r3; 776 c[colam + am3 + i] += r4; 777 } 778 b1 += bm4; 779 b2 += bm4; 780 b3 += bm4; 781 b4 += bm4; 782 } 783 for (;col<cn; col++){ /* over extra columns of C */ 784 for (i=0; i<am; i++) { /* over rows of C in those columns */ 785 r1 = 0.0; 786 n = a->i[i+1] - a->i[i]; 787 aj = a->j + a->i[i]; 788 aa = a->a + a->i[i]; 789 790 for (j=0; j<n; j++) { 791 r1 += (*aa++)*b1[*aj++]; 792 } 793 c[col*am + i] += r1; 794 } 795 b1 += bm; 796 } 797 } 798 ierr = PetscLogFlops(cn*2.0*a->nz);CHKERRQ(ierr); 799 ierr = MatRestoreArray(B,&b);CHKERRQ(ierr); 800 ierr = MatRestoreArray(C,&c);CHKERRQ(ierr); 801 PetscFunctionReturn(0); 802 } 803 804 #undef __FUNCT__ 805 #define __FUNCT__ "MatTransColoringApplySpToDen_SeqAIJ" 806 PetscErrorCode MatTransColoringApplySpToDen_SeqAIJ(MatTransposeColoring coloring,Mat B,Mat Btdense) 807 { 808 PetscErrorCode ierr; 809 Mat_SeqAIJ *a = (Mat_SeqAIJ*)B->data; 810 Mat_SeqDense *atdense = (Mat_SeqDense*)Btdense->data; 811 PetscInt m=Btdense->rmap->n,n=Btdense->cmap->n,j,k,l,col,anz,*atcol,brow; 812 MatScalar *atval,*bval=atdense->v; 813 PetscInt *columns=coloring->columns,*colorforcol=coloring->colorforcol; 814 815 PetscFunctionBegin; 816 ierr = PetscMemzero(atdense->v,(m*n)*sizeof(MatScalar));CHKERRQ(ierr); 817 for (k=0; k<coloring->ncolors; k++) { 818 for (l=0; l<coloring->ncolumns[k]; l++) { /* insert a row of B to a column of Btdense */ 819 col = *(columns + colorforcol[k] + l); 820 atcol = a->j + a->i[col]; 821 atval = a->a + a->i[col]; 822 anz = a->i[col+1] - a->i[col]; 823 for (j=0; j<anz; j++){ 824 brow = atcol[j]; 825 bval[k*m+brow] = atval[j]; 826 } 827 } 828 } 829 PetscFunctionReturn(0); 830 } 831 832 #undef __FUNCT__ 833 #define __FUNCT__ "MatTransColoringApplyDenToSp_SeqAIJ" 834 PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring,Mat Cden,Mat Csp) 835 { 836 PetscErrorCode ierr; 837 Mat_SeqAIJ *csp = (Mat_SeqAIJ*)Csp->data; 838 PetscInt k,l,row,m; 839 PetscScalar *ca_den,*cp_den,*ca=csp->a; 840 PetscInt *rows=matcoloring->rows,*spidx=matcoloring->columnsforspidx,*colorforrow=matcoloring->colorforrow,idx; 841 842 PetscFunctionBegin; 843 ierr = MatGetLocalSize(Csp,&m,PETSC_NULL);CHKERRQ(ierr); 844 ierr = MatGetArray(Cden,&ca_den);CHKERRQ(ierr); 845 cp_den = ca_den; 846 for (k=0; k<matcoloring->ncolors; k++) { 847 for (l=0; l<matcoloring->nrows[k]; l++){ 848 row = *(rows + colorforrow[k] + l); 849 idx = *(spidx + colorforrow[k] + l); 850 ca[idx] = cp_den[row]; 851 } 852 cp_den += m; 853 } 854 ierr = MatRestoreArray(Cden,&ca_den);CHKERRQ(ierr); 855 PetscFunctionReturn(0); 856 } 857 858 /* 859 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 860 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 861 spidx[], index of a->j, to be used for setting 'columnsforspidx' in MatTransposeColoringCreate_SeqAIJ(). 862 */ 863 #undef __FUNCT__ 864 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ_Color" 865 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 866 { 867 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 868 PetscErrorCode ierr; 869 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 870 PetscInt nz = a->i[m],row,*jj,mr,col; 871 PetscInt *cspidx; 872 873 PetscFunctionBegin; 874 *nn = n; 875 if (!ia) PetscFunctionReturn(0); 876 if (symmetric) { 877 SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"MatGetColumnIJ_SeqAIJ_Color() not supported for the case symmetric"); 878 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 879 } else { 880 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 881 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 882 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 883 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 884 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cspidx);CHKERRQ(ierr); 885 jj = a->j; 886 for (i=0; i<nz; i++) { 887 collengths[jj[i]]++; 888 } 889 cia[0] = oshift; 890 for (i=0; i<n; i++) { 891 cia[i+1] = cia[i] + collengths[i]; 892 } 893 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 894 jj = a->j; 895 for (row=0; row<m; row++) { 896 mr = a->i[row+1] - a->i[row]; 897 for (i=0; i<mr; i++) { 898 col = *jj++; 899 cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */ 900 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 901 } 902 } 903 ierr = PetscFree(collengths);CHKERRQ(ierr); 904 *ia = cia; *ja = cja; 905 *spidx = cspidx; 906 } 907 PetscFunctionReturn(0); 908 } 909 910 #undef __FUNCT__ 911 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ_Color" 912 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 913 { 914 PetscErrorCode ierr; 915 916 PetscFunctionBegin; 917 if (!ia) PetscFunctionReturn(0); 918 919 ierr = PetscFree(*ia);CHKERRQ(ierr); 920 ierr = PetscFree(*ja);CHKERRQ(ierr); 921 ierr = PetscFree(*spidx);CHKERRQ(ierr); 922 PetscFunctionReturn(0); 923 } 924 925 #undef __FUNCT__ 926 #define __FUNCT__ "MatTransposeColoringCreate_SeqAIJ" 927 PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatTransposeColoring c) 928 { 929 PetscErrorCode ierr; 930 PetscInt i,n,nrows,N,j,k,m,*row_idx,*ci,*cj,ncols,col,cm; 931 const PetscInt *is; 932 PetscInt nis = iscoloring->n,*rowhit,bs = 1; 933 IS *isa; 934 PetscBool done; 935 PetscBool flg1,flg2; 936 Mat_SeqAIJ *csp = (Mat_SeqAIJ*)mat->data; 937 PetscInt *colorforrow,*rows,*rows_i,*columnsforspidx,*columnsforspidx_i,*idxhit,*spidx; 938 PetscInt *colorforcol,*columns,*columns_i; 939 PetscLogDouble t0,tf,etime_GetIS=0.0,etime_GetColumnIJ=0.0,etime_memzero=0.0,etime_loop=0.0,t00,tf0; 940 941 PetscFunctionBegin; 942 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 943 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 944 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 945 etime_GetIS += tf - t0; 946 947 /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */ 948 ierr = PetscTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);CHKERRQ(ierr); 949 ierr = PetscTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);CHKERRQ(ierr); 950 if (flg1 || flg2) { 951 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 952 } 953 954 N = mat->cmap->N/bs; 955 c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */ 956 c->N = mat->cmap->N/bs; 957 c->m = mat->rmap->N/bs; 958 c->rstart = 0; 959 960 c->ncolors = nis; 961 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 962 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 963 ierr = PetscMalloc2(csp->nz+1,PetscInt,&rows,csp->nz+1,PetscInt,&columnsforspidx);CHKERRQ(ierr); 964 ierr = PetscMalloc((nis+1)*sizeof(PetscInt),&colorforrow);CHKERRQ(ierr); 965 colorforrow[0] = 0; 966 rows_i = rows; 967 columnsforspidx_i = columnsforspidx; 968 969 ierr = PetscMalloc((nis+1)*sizeof(PetscInt),&colorforcol);CHKERRQ(ierr); 970 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&columns);CHKERRQ(ierr); 971 colorforcol[0] = 0; 972 columns_i = columns; 973 974 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 975 ierr = MatGetColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,&done);CHKERRQ(ierr); /* column-wise storage of mat */ 976 if (!done) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"MatGetColumnIJ() not supported for matrix type %s",((PetscObject)mat)->type_name); 977 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 978 etime_GetColumnIJ += tf - t0; 979 980 cm = c->m; 981 ierr = PetscMalloc((cm+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 982 ierr = PetscMalloc((cm+1)*sizeof(PetscInt),&idxhit);CHKERRQ(ierr); 983 ierr = PetscGetTime(&t0);CHKERRQ(ierr); 984 for (i=0; i<nis; i++) { 985 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 986 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 987 c->ncolumns[i] = n; 988 if (n) { 989 ierr = PetscMemcpy(columns_i,is,n*sizeof(PetscInt));CHKERRQ(ierr); 990 } 991 colorforcol[i+1] = colorforcol[i] + n; 992 columns_i += n; 993 994 /* fast, crude version requires O(N*N) work */ 995 ierr = PetscGetTime(&t00);CHKERRQ(ierr); 996 ierr = PetscMemzero(rowhit,cm*sizeof(PetscInt));CHKERRQ(ierr); 997 ierr = PetscGetTime(&tf0);CHKERRQ(ierr); 998 etime_memzero += tf0 - t00; 999 1000 /* loop over columns*/ 1001 for (j=0; j<n; j++) { 1002 col = is[j]; 1003 row_idx = cj + ci[col]; 1004 m = ci[col+1] - ci[col]; 1005 /* loop over columns marking them in rowhit */ 1006 for (k=0; k<m; k++) { 1007 idxhit[*row_idx] = spidx[ci[col] + k]; 1008 rowhit[*row_idx++] = col + 1; 1009 } 1010 } 1011 /* count the number of hits */ 1012 nrows = 0; 1013 for (j=0; j<cm; j++) { 1014 if (rowhit[j]) nrows++; 1015 } 1016 c->nrows[i] = nrows; 1017 colorforrow[i+1] = colorforrow[i] + nrows; 1018 1019 nrows = 0; 1020 for (j=0; j<cm; j++) { 1021 if (rowhit[j]) { 1022 rows_i[nrows] = j; 1023 columnsforspidx_i[nrows] = idxhit[j]; 1024 nrows++; 1025 } 1026 } 1027 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 1028 rows_i += nrows; columnsforspidx_i += nrows; 1029 } 1030 ierr = PetscGetTime(&tf);CHKERRQ(ierr); 1031 etime_loop += tf - t0; 1032 1033 ierr = MatRestoreColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,&done);CHKERRQ(ierr); 1034 ierr = PetscFree(rowhit);CHKERRQ(ierr); 1035 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 1036 #if defined(PETSC_USE_DEBUG) 1037 if (csp->nz != colorforrow[nis]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"csp->nz %d != colorforrow[nis] %d",csp->nz,colorforrow[nis]); 1038 #endif 1039 1040 c->colorforrow = colorforrow; 1041 c->rows = rows; 1042 c->columnsforspidx = columnsforspidx; 1043 c->colorforcol = colorforcol; 1044 c->columns = columns; 1045 printf(" Time for ColoringCreate: IS %g, ColumnIJ: %g, memzero %g, loop %g\n",etime_GetIS,etime_GetColumnIJ,etime_memzero,etime_loop); 1046 ierr = PetscFree(idxhit); 1047 PetscFunctionReturn(0); 1048 } 1049