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