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