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