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(PETSC_SUCCESS); 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(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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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]) caj[k] += (*pA) * apa[apj[nextap++]]; 258 } 259 PetscCall(PetscLogFlops(2.0 * apnzj)); 260 pA++; 261 } 262 263 /* Zero the current row info for A*P */ 264 for (j = 0; j < apnzj; j++) { 265 apa[apj[j]] = 0.; 266 apjdense[apj[j]] = 0; 267 } 268 } 269 270 /* Assemble the final matrix and clean up */ 271 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 272 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 273 274 PetscCall(PetscFree2(apa, apjdense)); 275 PetscCall(PetscFree(apj)); 276 PetscFunctionReturn(PETSC_SUCCESS); 277 } 278 279 PetscErrorCode MatPtAPNumeric_SeqAIJ_SeqAIJ(Mat A, Mat P, Mat C) 280 { 281 Mat_MatTransMatMult *atb; 282 283 PetscFunctionBegin; 284 MatCheckProduct(C, 3); 285 atb = (Mat_MatTransMatMult *)C->product->data; 286 PetscCheck(atb, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing data structure"); 287 PetscCall(MatTranspose(P, MAT_REUSE_MATRIX, &atb->At)); 288 PetscCheck(C->ops->matmultnumeric, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing numeric operation"); 289 /* when using rap, MatMatMatMultSymbolic used a different data */ 290 if (atb->data) C->product->data = atb->data; 291 PetscCall((*C->ops->matmatmultnumeric)(atb->At, A, P, C)); 292 C->product->data = atb; 293 PetscFunctionReturn(PETSC_SUCCESS); 294 } 295