1 2 /* 3 Defines projective product routines where A is a SeqAIJ matrix 4 C = R * A * R^T 5 */ 6 7 #include <../src/mat/impls/aij/seq/aij.h> 8 #include <../src/mat/utils/freespace.h> 9 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatDestroy_SeqAIJ_RARt" 13 PetscErrorCode MatDestroy_SeqAIJ_RARt(Mat A) 14 { 15 PetscErrorCode ierr; 16 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 17 Mat_RARt *rart = a->rart; 18 19 PetscFunctionBegin; 20 ierr = MatTransposeColoringDestroy(&rart->matcoloring);CHKERRQ(ierr); 21 ierr = MatDestroy(&rart->Rt);CHKERRQ(ierr); 22 ierr = MatDestroy(&rart->RARt);CHKERRQ(ierr); 23 ierr = MatDestroy(&rart->ARt);CHKERRQ(ierr); 24 ierr = PetscFree(rart->work);CHKERRQ(ierr); 25 26 A->ops->destroy = rart->destroy; 27 if (A->ops->destroy) { 28 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 29 } 30 ierr = PetscFree(rart);CHKERRQ(ierr); 31 PetscFunctionReturn(0); 32 } 33 34 #undef __FUNCT__ 35 #define __FUNCT__ "MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart" 36 PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,PetscReal fill,Mat *C) 37 { 38 PetscErrorCode ierr; 39 Mat P; 40 PetscInt *rti,*rtj; 41 Mat_RARt *rart; 42 MatTransposeColoring matcoloring; 43 ISColoring iscoloring; 44 Mat Rt_dense,RARt_dense; 45 Mat_SeqAIJ *c; 46 47 PetscFunctionBegin; 48 /* create symbolic P=Rt */ 49 ierr = MatGetSymbolicTranspose_SeqAIJ(R,&rti,&rtj);CHKERRQ(ierr); 50 ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,R->cmap->n,R->rmap->n,rti,rtj,NULL,&P);CHKERRQ(ierr); 51 52 /* get symbolic C=Pt*A*P */ 53 ierr = MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A,P,fill,C);CHKERRQ(ierr); 54 (*C)->rmap->bs = R->rmap->bs; 55 (*C)->cmap->bs = R->rmap->bs; 56 (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart; 57 58 /* create a supporting struct */ 59 ierr = PetscNew(Mat_RARt,&rart);CHKERRQ(ierr); 60 c = (Mat_SeqAIJ*)(*C)->data; 61 c->rart = rart; 62 63 /* ------ Use coloring ---------- */ 64 /* inode causes memory problem, don't know why */ 65 if (c->inode.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MAT_USE_INODES is not supported. Use '-mat_no_inode'"); 66 67 /* Create MatTransposeColoring from symbolic C=R*A*R^T */ 68 ierr = MatGetColoring(*C,MATCOLORINGLF,&iscoloring);CHKERRQ(ierr); 69 ierr = MatTransposeColoringCreate(*C,iscoloring,&matcoloring);CHKERRQ(ierr); 70 71 rart->matcoloring = matcoloring; 72 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 73 74 /* Create Rt_dense */ 75 ierr = MatCreate(PETSC_COMM_SELF,&Rt_dense);CHKERRQ(ierr); 76 ierr = MatSetSizes(Rt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);CHKERRQ(ierr); 77 ierr = MatSetType(Rt_dense,MATSEQDENSE);CHKERRQ(ierr); 78 ierr = MatSeqDenseSetPreallocation(Rt_dense,NULL);CHKERRQ(ierr); 79 80 Rt_dense->assembled = PETSC_TRUE; 81 rart->Rt = Rt_dense; 82 83 /* Create RARt_dense = R*A*Rt_dense */ 84 ierr = MatCreate(PETSC_COMM_SELF,&RARt_dense);CHKERRQ(ierr); 85 ierr = MatSetSizes(RARt_dense,(*C)->rmap->n,matcoloring->ncolors,(*C)->rmap->n,matcoloring->ncolors);CHKERRQ(ierr); 86 ierr = MatSetType(RARt_dense,MATSEQDENSE);CHKERRQ(ierr); 87 ierr = MatSeqDenseSetPreallocation(RARt_dense,NULL);CHKERRQ(ierr); 88 89 rart->RARt = RARt_dense; 90 91 /* Allocate work array to store columns of A*R^T used in MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense() */ 92 ierr = PetscMalloc(A->rmap->n*4*sizeof(PetscScalar),&rart->work);CHKERRQ(ierr); 93 94 rart->destroy = (*C)->ops->destroy; 95 (*C)->ops->destroy = MatDestroy_SeqAIJ_RARt; 96 97 /* clean up */ 98 ierr = MatRestoreSymbolicTranspose_SeqAIJ(R,&rti,&rtj);CHKERRQ(ierr); 99 ierr = MatDestroy(&P);CHKERRQ(ierr); 100 101 #if defined(PETSC_USE_INFO) 102 { 103 PetscReal density= (PetscReal)(c->nz)/(RARt_dense->rmap->n*RARt_dense->cmap->n); 104 ierr = PetscInfo(*C,"C=R*(A*Rt) via coloring C - use sparse-dense inner products\n");CHKERRQ(ierr); 105 ierr = PetscInfo6(*C,"RARt_den %D %D; Rt %D %D (RARt->nz %D)/(m*ncolors)=%g\n",RARt_dense->rmap->n,RARt_dense->cmap->n,R->cmap->n,R->rmap->n,c->nz,density);CHKERRQ(ierr); 106 } 107 #endif 108 PetscFunctionReturn(0); 109 } 110 111 /* 112 RAB = R * A * B, R and A in seqaij format, B in dense format; 113 */ 114 #undef __FUNCT__ 115 #define __FUNCT__ "MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense" 116 PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(Mat R,Mat A,Mat B,Mat RAB,PetscScalar *work) 117 { 118 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*r=(Mat_SeqAIJ*)R->data; 119 PetscErrorCode ierr; 120 PetscScalar *b,r1,r2,r3,r4,*b1,*b2,*b3,*b4; 121 MatScalar *aa,*ra; 122 PetscInt cn =B->cmap->n,bm=B->rmap->n,col,i,j,n,*ai=a->i,*aj,am=A->rmap->n; 123 PetscInt am2=2*am,am3=3*am,bm4=4*bm; 124 PetscScalar *d,*c,*c2,*c3,*c4; 125 PetscInt *rj,rm=R->rmap->n,dm=RAB->rmap->n,dn=RAB->cmap->n; 126 PetscInt rm2=2*rm,rm3=3*rm,colrm; 127 128 PetscFunctionBegin; 129 if (!dm || !dn) PetscFunctionReturn(0); 130 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); 131 if (am != R->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in R %D not equal rows in A %D\n",R->cmap->n,am); 132 if (R->rmap->n != RAB->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows in RAB %D not equal rows in R %D\n",RAB->rmap->n,R->rmap->n); 133 if (B->cmap->n != RAB->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in RAB %D not equal columns in B %D\n",RAB->cmap->n,B->cmap->n); 134 135 { /* 136 This approach is not as good as original ones (will be removed later), but it reveals that 137 AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/examples/tutorials/ex56.c 138 */ 139 PetscBool via_matmatmult=PETSC_FALSE; 140 ierr = PetscOptionsGetBool(NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);CHKERRQ(ierr); 141 if (via_matmatmult) { 142 Mat AB_den; 143 ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,&AB_den);CHKERRQ(ierr); 144 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);CHKERRQ(ierr); 145 ierr = MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);CHKERRQ(ierr); 146 ierr = MatDestroy(&AB_den);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 } 150 151 ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr); 152 ierr = MatDenseGetArray(RAB,&d);CHKERRQ(ierr); 153 b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; 154 c = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am; 155 for (col=0; col<cn-4; col += 4) { /* over columns of C */ 156 for (i=0; i<am; i++) { /* over rows of A in those columns */ 157 r1 = r2 = r3 = r4 = 0.0; 158 n = ai[i+1] - ai[i]; 159 aj = a->j + ai[i]; 160 aa = a->a + ai[i]; 161 for (j=0; j<n; j++) { 162 r1 += (*aa)*b1[*aj]; 163 r2 += (*aa)*b2[*aj]; 164 r3 += (*aa)*b3[*aj]; 165 r4 += (*aa++)*b4[*aj++]; 166 } 167 c[i] = r1; 168 c[am + i] = r2; 169 c[am2 + i] = r3; 170 c[am3 + i] = r4; 171 } 172 b1 += bm4; 173 b2 += bm4; 174 b3 += bm4; 175 b4 += bm4; 176 177 /* RAB[:,col] = R*C[:,col] */ 178 colrm = col*rm; 179 for (i=0; i<rm; i++) { /* over rows of R in those columns */ 180 r1 = r2 = r3 = r4 = 0.0; 181 n = r->i[i+1] - r->i[i]; 182 rj = r->j + r->i[i]; 183 ra = r->a + r->i[i]; 184 for (j=0; j<n; j++) { 185 r1 += (*ra)*c[*rj]; 186 r2 += (*ra)*c2[*rj]; 187 r3 += (*ra)*c3[*rj]; 188 r4 += (*ra++)*c4[*rj++]; 189 } 190 d[colrm + i] = r1; 191 d[colrm + rm + i] = r2; 192 d[colrm + rm2 + i] = r3; 193 d[colrm + rm3 + i] = r4; 194 } 195 } 196 for (; col<cn; col++) { /* over extra columns of C */ 197 for (i=0; i<am; i++) { /* over rows of A in those columns */ 198 r1 = 0.0; 199 n = a->i[i+1] - a->i[i]; 200 aj = a->j + a->i[i]; 201 aa = a->a + a->i[i]; 202 for (j=0; j<n; j++) { 203 r1 += (*aa++)*b1[*aj++]; 204 } 205 c[i] = r1; 206 } 207 b1 += bm; 208 209 for (i=0; i<rm; i++) { /* over rows of R in those columns */ 210 r1 = 0.0; 211 n = r->i[i+1] - r->i[i]; 212 rj = r->j + r->i[i]; 213 ra = r->a + r->i[i]; 214 for (j=0; j<n; j++) { 215 r1 += (*ra++)*c[*rj++]; 216 } 217 d[col*rm + i] = r1; 218 } 219 } 220 ierr = PetscLogFlops(cn*2.0*(a->nz + r->nz));CHKERRQ(ierr); 221 222 ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr); 223 ierr = MatDenseRestoreArray(RAB,&d);CHKERRQ(ierr); 224 ierr = MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 225 ierr = MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 226 PetscFunctionReturn(0); 227 } 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart" 231 PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C) 232 { 233 PetscErrorCode ierr; 234 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 235 Mat_RARt *rart=c->rart; 236 MatTransposeColoring matcoloring; 237 Mat Rt,RARt; 238 239 PetscFunctionBegin; 240 /* Get dense Rt by Apply MatTransposeColoring to R */ 241 matcoloring = rart->matcoloring; 242 Rt = rart->Rt; 243 ierr = MatTransColoringApplySpToDen(matcoloring,R,Rt);CHKERRQ(ierr); 244 245 /* Get dense RARt = R*A*Rt -- dominates! */ 246 RARt = rart->RARt; 247 ierr = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);CHKERRQ(ierr); 248 249 /* Recover C from C_dense */ 250 ierr = MatTransColoringApplyDenToSp(matcoloring,RARt,C);CHKERRQ(ierr); 251 PetscFunctionReturn(0); 252 } 253 254 #undef __FUNCT__ 255 #define __FUNCT__ "MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult" 256 PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat *C) 257 { 258 PetscErrorCode ierr; 259 Mat ARt,RARt; 260 Mat_SeqAIJ *c; 261 Mat_RARt *rart; 262 263 PetscFunctionBegin; 264 /* must use '-mat_no_inode' with '-matmattransmult_color 1' - do not knwo why? */ 265 ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,R,fill,&ARt);CHKERRQ(ierr); 266 ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,&RARt);CHKERRQ(ierr); 267 *C = RARt; 268 RARt->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult; 269 270 ierr = PetscNew(Mat_RARt,&rart);CHKERRQ(ierr); 271 c = (Mat_SeqAIJ*)(*C)->data; 272 c->rart = rart; 273 rart->ARt = ARt; 274 rart->destroy = RARt->ops->destroy; 275 RARt->ops->destroy = MatDestroy_SeqAIJ_RARt; 276 #if defined(PETSC_USE_INFO) 277 ierr = PetscInfo(*C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");CHKERRQ(ierr); 278 #endif 279 PetscFunctionReturn(0); 280 } 281 282 #undef __FUNCT__ 283 #define __FUNCT__ "MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult" 284 PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C) 285 { 286 PetscErrorCode ierr; 287 Mat_SeqAIJ *c=(Mat_SeqAIJ*)C->data; 288 Mat_RARt *rart=c->rart; 289 Mat ARt=rart->ARt; 290 291 PetscFunctionBegin; 292 ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,ARt);CHKERRQ(ierr); /* dominate! */ 293 ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ(R,ARt,C);CHKERRQ(ierr); 294 PetscFunctionReturn(0); 295 } 296 297 #undef __FUNCT__ 298 #define __FUNCT__ "MatRARtSymbolic_SeqAIJ_SeqAIJ" 299 PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat R,PetscReal fill,Mat *C) 300 { 301 PetscErrorCode ierr; 302 Mat Rt; 303 Mat_SeqAIJ *c; 304 Mat_RARt *rart; 305 306 PetscFunctionBegin; 307 ierr = MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);CHKERRQ(ierr); 308 ierr = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);CHKERRQ(ierr); 309 310 ierr = PetscNew(Mat_RARt,&rart);CHKERRQ(ierr); 311 rart->Rt = Rt; 312 c = (Mat_SeqAIJ*)(*C)->data; 313 c->rart = rart; 314 rart->destroy = (*C)->ops->destroy; 315 (*C)->ops->destroy = MatDestroy_SeqAIJ_RARt; 316 (*C)->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ; 317 #if defined(PETSC_USE_INFO) 318 ierr = PetscInfo(*C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");CHKERRQ(ierr); 319 #endif 320 PetscFunctionReturn(0); 321 } 322 323 #undef __FUNCT__ 324 #define __FUNCT__ "MatRARtNumeric_SeqAIJ_SeqAIJ" 325 PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C) 326 { 327 PetscErrorCode ierr; 328 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 329 Mat_RARt *rart = c->rart; 330 Mat Rt = rart->Rt; 331 332 PetscFunctionBegin; 333 ierr = MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&Rt);CHKERRQ(ierr); 334 ierr = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,C);CHKERRQ(ierr); 335 PetscFunctionReturn(0); 336 } 337 338 #undef __FUNCT__ 339 #define __FUNCT__ "MatRARt_SeqAIJ_SeqAIJ" 340 PetscErrorCode MatRARt_SeqAIJ_SeqAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C) 341 { 342 PetscErrorCode ierr; 343 const char *algTypes[3] = {"matmatmatmult","matmattransposemult","coloring_rart"}; 344 PetscInt alg=0; /* set default algorithm */ 345 346 PetscFunctionBegin; 347 if (scall == MAT_INITIAL_MATRIX) { 348 ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr); 349 ierr = PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,3,algTypes[0],&alg,NULL);CHKERRQ(ierr); 350 ierr = PetscOptionsEnd();CHKERRQ(ierr); 351 352 ierr = PetscLogEventBegin(MAT_RARtSymbolic,A,R,0,0);CHKERRQ(ierr); 353 switch (alg) { 354 case 1: 355 /* via matmattransposemult: ARt=A*R^T, C=R*ARt - matrix coloring can be applied to A*R^T */ 356 ierr = MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,C);CHKERRQ(ierr); 357 break; 358 case 2: 359 /* via coloring_rart: apply coloring C = R*A*R^T */ 360 ierr = MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,C);CHKERRQ(ierr); 361 break; 362 default: 363 /* via matmatmatmult: Rt=R^T, C=R*A*Rt - avoid inefficient sparse inner products */ 364 ierr = MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,C);CHKERRQ(ierr); 365 break; 366 } 367 ierr = PetscLogEventEnd(MAT_RARtSymbolic,A,R,0,0);CHKERRQ(ierr); 368 } 369 370 ierr = PetscLogEventBegin(MAT_RARtNumeric,A,R,0,0);CHKERRQ(ierr); 371 ierr = (*(*C)->ops->rartnumeric)(A,R,*C);CHKERRQ(ierr); 372 ierr = PetscLogEventEnd(MAT_RARtNumeric,A,R,0,0);CHKERRQ(ierr); 373 PetscFunctionReturn(0); 374 } 375