1 2 /* 3 Defines projective product routines where A is a SeqAIJ matrix 4 C = P^T * A * P 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 <petsctime.h> 11 12 #if defined(PETSC_HAVE_HYPRE) 13 PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat,Mat,PetscReal,Mat); 14 #endif 15 16 PetscErrorCode MatProductSymbolic_PtAP_SeqAIJ_SeqAIJ(Mat C) 17 { 18 Mat_Product *product = C->product; 19 Mat A=product->A,P=product->B; 20 MatProductAlgorithm alg=product->alg; 21 PetscReal fill=product->fill; 22 PetscBool flg; 23 Mat Pt; 24 25 PetscFunctionBegin; 26 /* "scalable" */ 27 PetscCall(PetscStrcmp(alg,"scalable",&flg)); 28 if (flg) { 29 PetscCall(MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A,P,fill,C)); 30 C->ops->productnumeric = MatProductNumeric_PtAP; 31 PetscFunctionReturn(0); 32 } 33 34 /* "rap" */ 35 PetscCall(PetscStrcmp(alg,"rap",&flg)); 36 if (flg) { 37 Mat_MatTransMatMult *atb; 38 39 PetscCall(PetscNew(&atb)); 40 PetscCall(MatTranspose_SeqAIJ(P,MAT_INITIAL_MATRIX,&Pt)); 41 PetscCall(MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Pt,A,P,fill,C)); 42 43 atb->At = Pt; 44 atb->data = C->product->data; 45 atb->destroy = C->product->destroy; 46 C->product->data = atb; 47 C->product->destroy = MatDestroy_SeqAIJ_MatTransMatMult; 48 C->ops->ptapnumeric = MatPtAPNumeric_SeqAIJ_SeqAIJ; 49 C->ops->productnumeric = MatProductNumeric_PtAP; 50 PetscFunctionReturn(0); 51 } 52 53 /* hypre */ 54 #if defined(PETSC_HAVE_HYPRE) 55 PetscCall(PetscStrcmp(alg,"hypre",&flg)); 56 if (flg) { 57 PetscCall(MatPtAPSymbolic_AIJ_AIJ_wHYPRE(A,P,fill,C)); 58 PetscFunctionReturn(0); 59 } 60 #endif 61 62 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatProductType is not supported"); 63 } 64 65 PetscErrorCode MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(Mat A,Mat P,PetscReal fill,Mat C) 66 { 67 PetscFreeSpaceList free_space=NULL,current_space=NULL; 68 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*p = (Mat_SeqAIJ*)P->data,*c; 69 PetscInt *pti,*ptj,*ptJ,*ai=a->i,*aj=a->j,*ajj,*pi=p->i,*pj=p->j,*pjj; 70 PetscInt *ci,*cj,*ptadenserow,*ptasparserow,*ptaj,nspacedouble=0; 71 PetscInt an=A->cmap->N,am=A->rmap->N,pn=P->cmap->N,pm=P->rmap->N; 72 PetscInt i,j,k,ptnzi,arow,anzj,ptanzi,prow,pnzj,cnzi,nlnk,*lnk; 73 MatScalar *ca; 74 PetscBT lnkbt; 75 PetscReal afill; 76 77 PetscFunctionBegin; 78 /* Get ij structure of P^T */ 79 PetscCall(MatGetSymbolicTranspose_SeqAIJ(P,&pti,&ptj)); 80 ptJ = ptj; 81 82 /* Allocate ci array, arrays for fill computation and */ 83 /* free space for accumulating nonzero column info */ 84 PetscCall(PetscMalloc1(pn+1,&ci)); 85 ci[0] = 0; 86 87 PetscCall(PetscCalloc1(2*an+1,&ptadenserow)); 88 ptasparserow = ptadenserow + an; 89 90 /* create and initialize a linked list */ 91 nlnk = pn+1; 92 PetscCall(PetscLLCreate(pn,pn,nlnk,lnk,lnkbt)); 93 94 /* Set initial free space to be fill*(nnz(A)+ nnz(P)) */ 95 PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],pi[pm])),&free_space)); 96 current_space = free_space; 97 98 /* Determine symbolic info for each row of C: */ 99 for (i=0; i<pn; i++) { 100 ptnzi = pti[i+1] - pti[i]; 101 ptanzi = 0; 102 /* Determine symbolic row of PtA: */ 103 for (j=0; j<ptnzi; j++) { 104 arow = *ptJ++; 105 anzj = ai[arow+1] - ai[arow]; 106 ajj = aj + ai[arow]; 107 for (k=0; k<anzj; k++) { 108 if (!ptadenserow[ajj[k]]) { 109 ptadenserow[ajj[k]] = -1; 110 ptasparserow[ptanzi++] = ajj[k]; 111 } 112 } 113 } 114 /* Using symbolic info for row of PtA, determine symbolic info for row of C: */ 115 ptaj = ptasparserow; 116 cnzi = 0; 117 for (j=0; j<ptanzi; j++) { 118 prow = *ptaj++; 119 pnzj = pi[prow+1] - pi[prow]; 120 pjj = pj + pi[prow]; 121 /* add non-zero cols of P into the sorted linked list lnk */ 122 PetscCall(PetscLLAddSorted(pnzj,pjj,pn,&nlnk,lnk,lnkbt)); 123 cnzi += nlnk; 124 } 125 126 /* If free space is not available, make more free space */ 127 /* Double the amount of total space in the list */ 128 if (current_space->local_remaining<cnzi) { 129 PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),¤t_space)); 130 nspacedouble++; 131 } 132 133 /* Copy data into free space, and zero out denserows */ 134 PetscCall(PetscLLClean(pn,pn,cnzi,lnk,current_space->array,lnkbt)); 135 136 current_space->array += cnzi; 137 current_space->local_used += cnzi; 138 current_space->local_remaining -= cnzi; 139 140 for (j=0; j<ptanzi; j++) ptadenserow[ptasparserow[j]] = 0; 141 142 /* Aside: Perhaps we should save the pta info for the numerical factorization. */ 143 /* For now, we will recompute what is needed. */ 144 ci[i+1] = ci[i] + cnzi; 145 } 146 /* nnz is now stored in ci[ptm], column indices are in the list of free space */ 147 /* Allocate space for cj, initialize cj, and */ 148 /* destroy list of free space and other temporary array(s) */ 149 PetscCall(PetscMalloc1(ci[pn]+1,&cj)); 150 PetscCall(PetscFreeSpaceContiguous(&free_space,cj)); 151 PetscCall(PetscFree(ptadenserow)); 152 PetscCall(PetscLLDestroy(lnk,lnkbt)); 153 154 PetscCall(PetscCalloc1(ci[pn]+1,&ca)); 155 156 /* put together the new matrix */ 157 PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A),pn,pn,ci,cj,ca,((PetscObject)A)->type_name,C)); 158 PetscCall(MatSetBlockSizes(C,PetscAbs(P->cmap->bs),PetscAbs(P->cmap->bs))); 159 160 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ 161 /* Since these are PETSc arrays, change flags to free them as necessary. */ 162 c = (Mat_SeqAIJ*)((C)->data); 163 c->free_a = PETSC_TRUE; 164 c->free_ij = PETSC_TRUE; 165 c->nonew = 0; 166 167 C->ops->ptapnumeric = MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy; 168 169 /* set MatInfo */ 170 afill = (PetscReal)ci[pn]/(ai[am]+pi[pm] + 1.e-5); 171 if (afill < 1.0) afill = 1.0; 172 C->info.mallocs = nspacedouble; 173 C->info.fill_ratio_given = fill; 174 C->info.fill_ratio_needed = afill; 175 176 /* Clean up. */ 177 PetscCall(MatRestoreSymbolicTranspose_SeqAIJ(P,&pti,&ptj)); 178 #if defined(PETSC_USE_INFO) 179 if (ci[pn] != 0) { 180 PetscCall(PetscInfo(C,"Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n",nspacedouble,(double)fill,(double)afill)); 181 PetscCall(PetscInfo(C,"Use MatPtAP(A,P,MatReuse,%g,&C) for best performance.\n",(double)afill)); 182 } else { 183 PetscCall(PetscInfo(C,"Empty matrix product\n")); 184 } 185 #endif 186 PetscFunctionReturn(0); 187 } 188 189 PetscErrorCode MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy(Mat A,Mat P,Mat C) 190 { 191 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 192 Mat_SeqAIJ *p = (Mat_SeqAIJ*) P->data; 193 Mat_SeqAIJ *c = (Mat_SeqAIJ*) C->data; 194 PetscInt *ai=a->i,*aj=a->j,*apj,*apjdense,*pi=p->i,*pj=p->j,*pJ=p->j,*pjj; 195 PetscInt *ci=c->i,*cj=c->j,*cjj; 196 PetscInt am =A->rmap->N,cn=C->cmap->N,cm=C->rmap->N; 197 PetscInt i,j,k,anzi,pnzi,apnzj,nextap,pnzj,prow,crow; 198 MatScalar *aa,*apa,*pa,*pA,*paj,*ca,*caj; 199 200 PetscFunctionBegin; 201 /* Allocate temporary array for storage of one row of A*P (cn: non-scalable) */ 202 PetscCall(PetscCalloc2(cn,&apa,cn,&apjdense)); 203 PetscCall(PetscMalloc1(cn,&apj)); 204 /* trigger CPU copies if needed and flag CPU mask for C */ 205 #if defined(PETSC_HAVE_DEVICE) 206 { 207 const PetscScalar *dummy; 208 PetscCall(MatSeqAIJGetArrayRead(A,&dummy)); 209 PetscCall(MatSeqAIJRestoreArrayRead(A,&dummy)); 210 PetscCall(MatSeqAIJGetArrayRead(P,&dummy)); 211 PetscCall(MatSeqAIJRestoreArrayRead(P,&dummy)); 212 if (C->offloadmask != PETSC_OFFLOAD_UNALLOCATED) C->offloadmask = PETSC_OFFLOAD_CPU; 213 } 214 #endif 215 aa = a->a; 216 pa = p->a; 217 pA = p->a; 218 ca = c->a; 219 220 /* Clear old values in C */ 221 PetscCall(PetscArrayzero(ca,ci[cm])); 222 223 for (i=0; i<am; i++) { 224 /* Form sparse row of A*P */ 225 anzi = ai[i+1] - ai[i]; 226 apnzj = 0; 227 for (j=0; j<anzi; j++) { 228 prow = *aj++; 229 pnzj = pi[prow+1] - pi[prow]; 230 pjj = pj + pi[prow]; 231 paj = pa + pi[prow]; 232 for (k=0; k<pnzj; k++) { 233 if (!apjdense[pjj[k]]) { 234 apjdense[pjj[k]] = -1; 235 apj[apnzj++] = pjj[k]; 236 } 237 apa[pjj[k]] += (*aa)*paj[k]; 238 } 239 PetscCall(PetscLogFlops(2.0*pnzj)); 240 aa++; 241 } 242 243 /* Sort the j index array for quick sparse axpy. */ 244 /* Note: a array does not need sorting as it is in dense storage locations. */ 245 PetscCall(PetscSortInt(apnzj,apj)); 246 247 /* Compute P^T*A*P using outer product (P^T)[:,j]*(A*P)[j,:]. */ 248 pnzi = pi[i+1] - pi[i]; 249 for (j=0; j<pnzi; j++) { 250 nextap = 0; 251 crow = *pJ++; 252 cjj = cj + ci[crow]; 253 caj = ca + ci[crow]; 254 /* Perform sparse axpy operation. Note cjj includes apj. */ 255 for (k=0; nextap<apnzj; k++) { 256 PetscAssert(k < ci[crow+1] - ci[crow],PETSC_COMM_SELF,PETSC_ERR_PLIB,"k too large k %" PetscInt_FMT ", crow %" PetscInt_FMT,k,crow); 257 if (cjj[k]==apj[nextap]) { 258 caj[k] += (*pA)*apa[apj[nextap++]]; 259 } 260 } 261 PetscCall(PetscLogFlops(2.0*apnzj)); 262 pA++; 263 } 264 265 /* Zero the current row info for A*P */ 266 for (j=0; j<apnzj; j++) { 267 apa[apj[j]] = 0.; 268 apjdense[apj[j]] = 0; 269 } 270 } 271 272 /* Assemble the final matrix and clean up */ 273 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 274 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 275 276 PetscCall(PetscFree2(apa,apjdense)); 277 PetscCall(PetscFree(apj)); 278 PetscFunctionReturn(0); 279 } 280 281 PetscErrorCode MatPtAPNumeric_SeqAIJ_SeqAIJ(Mat A,Mat P,Mat C) 282 { 283 Mat_MatTransMatMult *atb; 284 285 PetscFunctionBegin; 286 MatCheckProduct(C,3); 287 atb = (Mat_MatTransMatMult*)C->product->data; 288 PetscCheck(atb,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing data structure"); 289 PetscCall(MatTranspose_SeqAIJ(P,MAT_REUSE_MATRIX,&atb->At)); 290 PetscCheck(C->ops->matmultnumeric,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing numeric operation"); 291 /* when using rap, MatMatMatMultSymbolic used a different data */ 292 if (atb->data) C->product->data = atb->data; 293 PetscCall((*C->ops->matmatmultnumeric)(atb->At,A,P,C)); 294 C->product->data = atb; 295 PetscFunctionReturn(0); 296 } 297