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