/* Defines matrix-matrix product routines for pairs of SeqAIJ matrices C = A * B */ #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ #include <../src/mat/utils/freespace.h> #include #include #include <../src/mat/impls/dense/seq/dense.h> PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C) { PetscFunctionBegin; if (C->ops->matmultnumeric) PetscCall((*C->ops->matmultnumeric)(A, B, C)); else PetscCall(MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(A, B, C)); PetscFunctionReturn(PETSC_SUCCESS); } /* Modified from MatCreateSeqAIJWithArrays() */ PETSC_INTERN PetscErrorCode MatSetSeqAIJWithArrays_private(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], MatType mtype, Mat mat) { PetscInt ii; Mat_SeqAIJ *aij; PetscBool isseqaij, ofree_a, ofree_ij; PetscFunctionBegin; PetscCheck(m <= 0 || !i[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0"); PetscCall(MatSetSizes(mat, m, n, m, n)); if (!mtype) { PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQAIJ, &isseqaij)); if (!isseqaij) PetscCall(MatSetType(mat, MATSEQAIJ)); } else { PetscCall(MatSetType(mat, mtype)); } aij = (Mat_SeqAIJ *)mat->data; ofree_a = aij->free_a; ofree_ij = aij->free_ij; /* changes the free flags */ PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat, MAT_SKIP_ALLOCATION, NULL)); PetscCall(PetscFree(aij->ilen)); PetscCall(PetscFree(aij->imax)); PetscCall(PetscMalloc1(m, &aij->imax)); PetscCall(PetscMalloc1(m, &aij->ilen)); for (ii = 0, aij->nonzerorowcnt = 0, aij->rmax = 0; ii < m; ii++) { const PetscInt rnz = i[ii + 1] - i[ii]; aij->nonzerorowcnt += !!rnz; aij->rmax = PetscMax(aij->rmax, rnz); aij->ilen[ii] = aij->imax[ii] = i[ii + 1] - i[ii]; } aij->maxnz = i[m]; aij->nz = i[m]; if (ofree_a) PetscCall(PetscShmgetDeallocateArray((void **)&aij->a)); if (ofree_ij) PetscCall(PetscShmgetDeallocateArray((void **)&aij->j)); if (ofree_ij) PetscCall(PetscShmgetDeallocateArray((void **)&aij->i)); aij->i = i; aij->j = j; aij->a = a; aij->nonew = -1; /* this indicates that inserting a new value in the matrix that generates a new nonzero is an error */ aij->free_a = PETSC_FALSE; aij->free_ij = PETSC_FALSE; PetscCall(MatCheckCompressedRow(mat, aij->nonzerorowcnt, &aij->compressedrow, aij->i, m, 0.6)); // Always build the diag info when i, j are set PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C) { Mat_Product *product = C->product; MatProductAlgorithm alg; PetscBool flg; PetscFunctionBegin; if (product) { alg = product->alg; } else { alg = "sorted"; } /* sorted */ PetscCall(PetscStrcmp(alg, "sorted", &flg)); if (flg) { PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(A, B, fill, C)); PetscFunctionReturn(PETSC_SUCCESS); } /* scalable */ PetscCall(PetscStrcmp(alg, "scalable", &flg)); if (flg) { PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A, B, fill, C)); PetscFunctionReturn(PETSC_SUCCESS); } /* scalable_fast */ PetscCall(PetscStrcmp(alg, "scalable_fast", &flg)); if (flg) { PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(A, B, fill, C)); PetscFunctionReturn(PETSC_SUCCESS); } /* heap */ PetscCall(PetscStrcmp(alg, "heap", &flg)); if (flg) { PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(A, B, fill, C)); PetscFunctionReturn(PETSC_SUCCESS); } /* btheap */ PetscCall(PetscStrcmp(alg, "btheap", &flg)); if (flg) { PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(A, B, fill, C)); PetscFunctionReturn(PETSC_SUCCESS); } /* llcondensed */ PetscCall(PetscStrcmp(alg, "llcondensed", &flg)); if (flg) { PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(A, B, fill, C)); PetscFunctionReturn(PETSC_SUCCESS); } /* rowmerge */ PetscCall(PetscStrcmp(alg, "rowmerge", &flg)); if (flg) { PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(A, B, fill, C)); PetscFunctionReturn(PETSC_SUCCESS); } #if defined(PETSC_HAVE_HYPRE) PetscCall(PetscStrcmp(alg, "hypre", &flg)); if (flg) { PetscCall(MatMatMultSymbolic_AIJ_AIJ_wHYPRE(A, B, fill, C)); PetscFunctionReturn(PETSC_SUCCESS); } #endif SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Mat Product Algorithm is not supported"); } PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(Mat A, Mat B, PetscReal fill, Mat C) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c; PetscInt *ai = a->i, *bi = b->i, *ci, *cj; PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N; PetscReal afill; PetscInt i, j, anzi, brow, bnzj, cnzi, *bj, *aj, *lnk, ndouble = 0, Crmax; PetscHMapI ta; PetscBT lnkbt; PetscFreeSpaceList free_space = NULL, current_space = NULL; PetscFunctionBegin; /* Get ci and cj */ /* Allocate ci array, arrays for fill computation and */ /* free space for accumulating nonzero column info */ PetscCall(PetscMalloc1(am + 2, &ci)); ci[0] = 0; /* create and initialize a linked list */ PetscCall(PetscHMapICreateWithSize(bn, &ta)); MatRowMergeMax_SeqAIJ(b, bm, ta); PetscCall(PetscHMapIGetSize(ta, &Crmax)); PetscCall(PetscHMapIDestroy(&ta)); PetscCall(PetscLLCondensedCreate(Crmax, bn, &lnk, &lnkbt)); /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space)); current_space = free_space; /* Determine ci and cj */ for (i = 0; i < am; i++) { anzi = ai[i + 1] - ai[i]; aj = a->j + ai[i]; for (j = 0; j < anzi; j++) { brow = aj[j]; bnzj = bi[brow + 1] - bi[brow]; bj = b->j + bi[brow]; /* add non-zero cols of B into the sorted linked list lnk */ PetscCall(PetscLLCondensedAddSorted(bnzj, bj, lnk, lnkbt)); } /* add possible missing diagonal entry */ if (C->force_diagonals) PetscCall(PetscLLCondensedAddSorted(1, &i, lnk, lnkbt)); cnzi = lnk[0]; /* If free space is not available, make more free space */ /* Double the amount of total space in the list */ if (current_space->local_remaining < cnzi) { PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), ¤t_space)); ndouble++; } /* Copy data into free space, then initialize lnk */ PetscCall(PetscLLCondensedClean(bn, cnzi, current_space->array, lnk, lnkbt)); current_space->array += cnzi; current_space->local_used += cnzi; current_space->local_remaining -= cnzi; ci[i + 1] = ci[i] + cnzi; } /* Column indices are in the list of free space */ /* Allocate space for cj, initialize cj, and */ /* destroy list of free space and other temporary array(s) */ PetscCall(PetscMalloc1(ci[am] + 1, &cj)); PetscCall(PetscFreeSpaceContiguous(&free_space, cj)); PetscCall(PetscLLCondensedDestroy(lnk, lnkbt)); /* put together the new symbolic matrix */ PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C)); PetscCall(MatSetBlockSizesFromMats(C, A, B)); /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ c = (Mat_SeqAIJ *)C->data; c->free_a = PETSC_FALSE; c->free_ij = PETSC_TRUE; c->nonew = 0; /* fast, needs non-scalable O(bn) array 'abdense' */ C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted; /* set MatInfo */ afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5; if (afill < 1.0) afill = 1.0; C->info.mallocs = ndouble; C->info.fill_ratio_given = fill; C->info.fill_ratio_needed = afill; #if defined(PETSC_USE_INFO) if (ci[am]) { PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill)); PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill)); } else { PetscCall(PetscInfo(C, "Empty matrix product\n")); } #endif PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(Mat A, Mat B, Mat C) { PetscLogDouble flops = 0.0; Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *bjj, *ci = c->i, *cj = c->j; PetscInt am = A->rmap->n, cm = C->rmap->n; PetscInt i, j, k, anzi, bnzi, cnzi, brow; PetscScalar *ca, valtmp; PetscScalar *ab_dense; PetscContainer cab_dense; const PetscScalar *aa, *ba, *baj; PetscFunctionBegin; PetscCall(MatSeqAIJGetArrayRead(A, &aa)); PetscCall(MatSeqAIJGetArrayRead(B, &ba)); if (!c->a) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */ PetscCall(PetscMalloc1(ci[cm] + 1, &ca)); c->a = ca; c->free_a = PETSC_TRUE; } else ca = c->a; /* TODO this should be done in the symbolic phase */ /* However, this function is so heavily used (sometimes in an hidden way through multnumeric function pointers that is hard to eradicate) */ PetscCall(PetscObjectQuery((PetscObject)C, "__PETSc__ab_dense", (PetscObject *)&cab_dense)); if (!cab_dense) { PetscCall(PetscMalloc1(B->cmap->N, &ab_dense)); PetscCall(PetscObjectContainerCompose((PetscObject)C, "__PETSc__ab_dense", ab_dense, PetscCtxDestroyDefault)); } else PetscCall(PetscContainerGetPointer(cab_dense, &ab_dense)); PetscCall(PetscArrayzero(ab_dense, B->cmap->N)); /* clean old values in C */ PetscCall(PetscArrayzero(ca, ci[cm])); /* Traverse A row-wise. */ /* Build the ith row in C by summing over nonzero columns in A, */ /* the rows of B corresponding to nonzeros of A. */ for (i = 0; i < am; i++) { anzi = ai[i + 1] - ai[i]; for (j = 0; j < anzi; j++) { brow = aj[j]; bnzi = bi[brow + 1] - bi[brow]; bjj = PetscSafePointerPlusOffset(bj, bi[brow]); baj = PetscSafePointerPlusOffset(ba, bi[brow]); /* perform dense axpy */ valtmp = aa[j]; for (k = 0; k < bnzi; k++) ab_dense[bjj[k]] += valtmp * baj[k]; flops += 2 * bnzi; } aj = PetscSafePointerPlusOffset(aj, anzi); aa = PetscSafePointerPlusOffset(aa, anzi); cnzi = ci[i + 1] - ci[i]; for (k = 0; k < cnzi; k++) { ca[k] += ab_dense[cj[k]]; ab_dense[cj[k]] = 0.0; /* zero ab_dense */ } flops += cnzi; cj = PetscSafePointerPlusOffset(cj, cnzi); ca += cnzi; } #if defined(PETSC_HAVE_DEVICE) if (C->offloadmask != PETSC_OFFLOAD_UNALLOCATED) C->offloadmask = PETSC_OFFLOAD_CPU; #endif PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); PetscCall(PetscLogFlops(flops)); PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(Mat A, Mat B, Mat C) { PetscLogDouble flops = 0.0; Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *bjj, *ci = c->i, *cj = c->j; PetscInt am = A->rmap->N, cm = C->rmap->N; PetscInt i, j, k, anzi, bnzi, cnzi, brow; PetscScalar *ca = c->a, valtmp; const PetscScalar *aa, *ba, *baj; PetscInt nextb; PetscFunctionBegin; PetscCall(MatSeqAIJGetArrayRead(A, &aa)); PetscCall(MatSeqAIJGetArrayRead(B, &ba)); if (!ca) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */ PetscCall(PetscMalloc1(ci[cm] + 1, &ca)); c->a = ca; c->free_a = PETSC_TRUE; } /* clean old values in C */ PetscCall(PetscArrayzero(ca, ci[cm])); /* Traverse A row-wise. */ /* Build the ith row in C by summing over nonzero columns in A, */ /* the rows of B corresponding to nonzeros of A. */ for (i = 0; i < am; i++) { anzi = ai[i + 1] - ai[i]; cnzi = ci[i + 1] - ci[i]; for (j = 0; j < anzi; j++) { brow = aj[j]; bnzi = bi[brow + 1] - bi[brow]; bjj = bj + bi[brow]; baj = ba + bi[brow]; /* perform sparse axpy */ valtmp = aa[j]; nextb = 0; for (k = 0; nextb < bnzi; k++) { if (cj[k] == bjj[nextb]) { /* ccol == bcol */ ca[k] += valtmp * baj[nextb++]; } } flops += 2 * bnzi; } aj += anzi; aa += anzi; cj += cnzi; ca += cnzi; } #if defined(PETSC_HAVE_DEVICE) if (C->offloadmask != PETSC_OFFLOAD_UNALLOCATED) C->offloadmask = PETSC_OFFLOAD_CPU; #endif PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); PetscCall(PetscLogFlops(flops)); PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(Mat A, Mat B, PetscReal fill, Mat C) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c; PetscInt *ai = a->i, *bi = b->i, *ci, *cj; PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N; MatScalar *ca; PetscReal afill; PetscInt i, j, anzi, brow, bnzj, cnzi, *bj, *aj, *lnk, ndouble = 0, Crmax; PetscHMapI ta; PetscFreeSpaceList free_space = NULL, current_space = NULL; PetscFunctionBegin; /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_fast() */ /* Allocate arrays for fill computation and free space for accumulating nonzero column */ PetscCall(PetscMalloc1(am + 2, &ci)); ci[0] = 0; /* create and initialize a linked list */ PetscCall(PetscHMapICreateWithSize(bn, &ta)); MatRowMergeMax_SeqAIJ(b, bm, ta); PetscCall(PetscHMapIGetSize(ta, &Crmax)); PetscCall(PetscHMapIDestroy(&ta)); PetscCall(PetscLLCondensedCreate_fast(Crmax, &lnk)); /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space)); current_space = free_space; /* Determine ci and cj */ for (i = 0; i < am; i++) { anzi = ai[i + 1] - ai[i]; aj = a->j + ai[i]; for (j = 0; j < anzi; j++) { brow = aj[j]; bnzj = bi[brow + 1] - bi[brow]; bj = b->j + bi[brow]; /* add non-zero cols of B into the sorted linked list lnk */ PetscCall(PetscLLCondensedAddSorted_fast(bnzj, bj, lnk)); } /* add possible missing diagonal entry */ if (C->force_diagonals) PetscCall(PetscLLCondensedAddSorted_fast(1, &i, lnk)); cnzi = lnk[1]; /* If free space is not available, make more free space */ /* Double the amount of total space in the list */ if (current_space->local_remaining < cnzi) { PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), ¤t_space)); ndouble++; } /* Copy data into free space, then initialize lnk */ PetscCall(PetscLLCondensedClean_fast(cnzi, current_space->array, lnk)); current_space->array += cnzi; current_space->local_used += cnzi; current_space->local_remaining -= cnzi; ci[i + 1] = ci[i] + cnzi; } /* Column indices are in the list of free space */ /* Allocate space for cj, initialize cj, and */ /* destroy list of free space and other temporary array(s) */ PetscCall(PetscMalloc1(ci[am] + 1, &cj)); PetscCall(PetscFreeSpaceContiguous(&free_space, cj)); PetscCall(PetscLLCondensedDestroy_fast(lnk)); /* Allocate space for ca */ PetscCall(PetscCalloc1(ci[am] + 1, &ca)); /* put together the new symbolic matrix */ PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, ca, ((PetscObject)A)->type_name, C)); PetscCall(MatSetBlockSizesFromMats(C, A, B)); /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ c = (Mat_SeqAIJ *)C->data; c->free_a = PETSC_TRUE; c->free_ij = PETSC_TRUE; c->nonew = 0; /* slower, less memory */ C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* set MatInfo */ afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5; if (afill < 1.0) afill = 1.0; C->info.mallocs = ndouble; C->info.fill_ratio_given = fill; C->info.fill_ratio_needed = afill; #if defined(PETSC_USE_INFO) if (ci[am]) { PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill)); PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill)); } else { PetscCall(PetscInfo(C, "Empty matrix product\n")); } #endif PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(Mat A, Mat B, PetscReal fill, Mat C) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c; PetscInt *ai = a->i, *bi = b->i, *ci, *cj; PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N; MatScalar *ca; PetscReal afill; PetscInt i, j, anzi, brow, bnzj, cnzi, *bj, *aj, *lnk, ndouble = 0, Crmax; PetscHMapI ta; PetscFreeSpaceList free_space = NULL, current_space = NULL; PetscFunctionBegin; /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_Scalalbe() */ /* Allocate arrays for fill computation and free space for accumulating nonzero column */ PetscCall(PetscMalloc1(am + 2, &ci)); ci[0] = 0; /* create and initialize a linked list */ PetscCall(PetscHMapICreateWithSize(bn, &ta)); MatRowMergeMax_SeqAIJ(b, bm, ta); PetscCall(PetscHMapIGetSize(ta, &Crmax)); PetscCall(PetscHMapIDestroy(&ta)); PetscCall(PetscLLCondensedCreate_Scalable(Crmax, &lnk)); /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space)); current_space = free_space; /* Determine ci and cj */ for (i = 0; i < am; i++) { anzi = ai[i + 1] - ai[i]; aj = a->j + ai[i]; for (j = 0; j < anzi; j++) { brow = aj[j]; bnzj = bi[brow + 1] - bi[brow]; bj = b->j + bi[brow]; /* add non-zero cols of B into the sorted linked list lnk */ PetscCall(PetscLLCondensedAddSorted_Scalable(bnzj, bj, lnk)); } /* add possible missing diagonal entry */ if (C->force_diagonals) PetscCall(PetscLLCondensedAddSorted_Scalable(1, &i, lnk)); cnzi = lnk[0]; /* If free space is not available, make more free space */ /* Double the amount of total space in the list */ if (current_space->local_remaining < cnzi) { PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), ¤t_space)); ndouble++; } /* Copy data into free space, then initialize lnk */ PetscCall(PetscLLCondensedClean_Scalable(cnzi, current_space->array, lnk)); current_space->array += cnzi; current_space->local_used += cnzi; current_space->local_remaining -= cnzi; ci[i + 1] = ci[i] + cnzi; } /* Column indices are in the list of free space */ /* Allocate space for cj, initialize cj, and */ /* destroy list of free space and other temporary array(s) */ PetscCall(PetscMalloc1(ci[am] + 1, &cj)); PetscCall(PetscFreeSpaceContiguous(&free_space, cj)); PetscCall(PetscLLCondensedDestroy_Scalable(lnk)); /* Allocate space for ca */ PetscCall(PetscCalloc1(ci[am] + 1, &ca)); /* put together the new symbolic matrix */ PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, ca, ((PetscObject)A)->type_name, C)); PetscCall(MatSetBlockSizesFromMats(C, A, B)); /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ c = (Mat_SeqAIJ *)C->data; c->free_a = PETSC_TRUE; c->free_ij = PETSC_TRUE; c->nonew = 0; /* slower, less memory */ C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* set MatInfo */ afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5; if (afill < 1.0) afill = 1.0; C->info.mallocs = ndouble; C->info.fill_ratio_given = fill; C->info.fill_ratio_needed = afill; #if defined(PETSC_USE_INFO) if (ci[am]) { PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill)); PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill)); } else { PetscCall(PetscInfo(C, "Empty matrix product\n")); } #endif PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(Mat A, Mat B, PetscReal fill, Mat C) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c; const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j; PetscInt *ci, *cj, *bb; PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N; PetscReal afill; PetscInt i, j, col, ndouble = 0; PetscFreeSpaceList free_space = NULL, current_space = NULL; PetscHeap h; PetscFunctionBegin; /* Get ci and cj - by merging sorted rows using a heap */ /* Allocate arrays for fill computation and free space for accumulating nonzero column */ PetscCall(PetscMalloc1(am + 2, &ci)); ci[0] = 0; /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space)); current_space = free_space; PetscCall(PetscHeapCreate(a->rmax, &h)); PetscCall(PetscMalloc1(a->rmax, &bb)); /* Determine ci and cj */ for (i = 0; i < am; i++) { const PetscInt anzi = ai[i + 1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */ const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */ ci[i + 1] = ci[i]; /* Populate the min heap */ for (j = 0; j < anzi; j++) { bb[j] = bi[acol[j]]; /* bb points at the start of the row */ if (bb[j] < bi[acol[j] + 1]) { /* Add if row is nonempty */ PetscCall(PetscHeapAdd(h, j, bj[bb[j]++])); } } /* Pick off the min element, adding it to free space */ PetscCall(PetscHeapPop(h, &j, &col)); while (j >= 0) { if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */ PetscCall(PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2, current_space->total_array_size), 16 << 20), ¤t_space)); ndouble++; } *(current_space->array++) = col; current_space->local_used++; current_space->local_remaining--; ci[i + 1]++; /* stash if anything else remains in this row of B */ if (bb[j] < bi[acol[j] + 1]) PetscCall(PetscHeapStash(h, j, bj[bb[j]++])); while (1) { /* pop and stash any other rows of B that also had an entry in this column */ PetscInt j2, col2; PetscCall(PetscHeapPeek(h, &j2, &col2)); if (col2 != col) break; PetscCall(PetscHeapPop(h, &j2, &col2)); if (bb[j2] < bi[acol[j2] + 1]) PetscCall(PetscHeapStash(h, j2, bj[bb[j2]++])); } /* Put any stashed elements back into the min heap */ PetscCall(PetscHeapUnstash(h)); PetscCall(PetscHeapPop(h, &j, &col)); } } PetscCall(PetscFree(bb)); PetscCall(PetscHeapDestroy(&h)); /* Column indices are in the list of free space */ /* Allocate space for cj, initialize cj, and */ /* destroy list of free space and other temporary array(s) */ PetscCall(PetscMalloc1(ci[am], &cj)); PetscCall(PetscFreeSpaceContiguous(&free_space, cj)); /* put together the new symbolic matrix */ PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C)); PetscCall(MatSetBlockSizesFromMats(C, A, B)); /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ c = (Mat_SeqAIJ *)C->data; c->free_a = PETSC_TRUE; c->free_ij = PETSC_TRUE; c->nonew = 0; C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted; /* set MatInfo */ afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5; if (afill < 1.0) afill = 1.0; C->info.mallocs = ndouble; C->info.fill_ratio_given = fill; C->info.fill_ratio_needed = afill; #if defined(PETSC_USE_INFO) if (ci[am]) { PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill)); PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill)); } else { PetscCall(PetscInfo(C, "Empty matrix product\n")); } #endif PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(Mat A, Mat B, PetscReal fill, Mat C) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c; const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j; PetscInt *ci, *cj, *bb; PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N; PetscReal afill; PetscInt i, j, col, ndouble = 0; PetscFreeSpaceList free_space = NULL, current_space = NULL; PetscHeap h; PetscBT bt; PetscFunctionBegin; /* Get ci and cj - using a heap for the sorted rows, but use BT so that each index is only added once */ /* Allocate arrays for fill computation and free space for accumulating nonzero column */ PetscCall(PetscMalloc1(am + 2, &ci)); ci[0] = 0; /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space)); current_space = free_space; PetscCall(PetscHeapCreate(a->rmax, &h)); PetscCall(PetscMalloc1(a->rmax, &bb)); PetscCall(PetscBTCreate(bn, &bt)); /* Determine ci and cj */ for (i = 0; i < am; i++) { const PetscInt anzi = ai[i + 1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */ const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */ const PetscInt *fptr = current_space->array; /* Save beginning of the row so we can clear the BT later */ ci[i + 1] = ci[i]; /* Populate the min heap */ for (j = 0; j < anzi; j++) { PetscInt brow = acol[j]; for (bb[j] = bi[brow]; bb[j] < bi[brow + 1]; bb[j]++) { PetscInt bcol = bj[bb[j]]; if (!PetscBTLookupSet(bt, bcol)) { /* new entry */ PetscCall(PetscHeapAdd(h, j, bcol)); bb[j]++; break; } } } /* Pick off the min element, adding it to free space */ PetscCall(PetscHeapPop(h, &j, &col)); while (j >= 0) { if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */ fptr = NULL; /* need PetscBTMemzero */ PetscCall(PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2, current_space->total_array_size), 16 << 20), ¤t_space)); ndouble++; } *(current_space->array++) = col; current_space->local_used++; current_space->local_remaining--; ci[i + 1]++; /* stash if anything else remains in this row of B */ for (; bb[j] < bi[acol[j] + 1]; bb[j]++) { PetscInt bcol = bj[bb[j]]; if (!PetscBTLookupSet(bt, bcol)) { /* new entry */ PetscCall(PetscHeapAdd(h, j, bcol)); bb[j]++; break; } } PetscCall(PetscHeapPop(h, &j, &col)); } if (fptr) { /* Clear the bits for this row */ for (; fptr < current_space->array; fptr++) PetscCall(PetscBTClear(bt, *fptr)); } else { /* We reallocated so we don't remember (easily) how to clear only the bits we changed */ PetscCall(PetscBTMemzero(bn, bt)); } } PetscCall(PetscFree(bb)); PetscCall(PetscHeapDestroy(&h)); PetscCall(PetscBTDestroy(&bt)); /* Column indices are in the list of free space */ /* Allocate space for cj, initialize cj, and */ /* destroy list of free space and other temporary array(s) */ PetscCall(PetscMalloc1(ci[am], &cj)); PetscCall(PetscFreeSpaceContiguous(&free_space, cj)); /* put together the new symbolic matrix */ PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C)); PetscCall(MatSetBlockSizesFromMats(C, A, B)); /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ c = (Mat_SeqAIJ *)C->data; c->free_a = PETSC_TRUE; c->free_ij = PETSC_TRUE; c->nonew = 0; C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted; /* set MatInfo */ afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5; if (afill < 1.0) afill = 1.0; C->info.mallocs = ndouble; C->info.fill_ratio_given = fill; C->info.fill_ratio_needed = afill; #if defined(PETSC_USE_INFO) if (ci[am]) { PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill)); PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill)); } else { PetscCall(PetscInfo(C, "Empty matrix product\n")); } #endif PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(Mat A, Mat B, PetscReal fill, Mat C) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c; const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j, *inputi, *inputj, *inputcol, *inputcol_L1; PetscInt *ci, *cj, *outputj, worki_L1[9], worki_L2[9]; PetscInt c_maxmem, a_maxrownnz = 0, a_rownnz; const PetscInt workcol[8] = {0, 1, 2, 3, 4, 5, 6, 7}; const PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N; const PetscInt *brow_ptr[8], *brow_end[8]; PetscInt window[8]; PetscInt window_min, old_window_min, ci_nnz, outputi_nnz = 0, L1_nrows, L2_nrows; PetscInt i, k, ndouble = 0, L1_rowsleft, rowsleft; PetscReal afill; PetscInt *workj_L1, *workj_L2, *workj_L3; PetscInt L1_nnz, L2_nnz; /* Step 1: Get upper bound on memory required for allocation. Because of the way virtual memory works, only the memory pages that are actually needed will be physically allocated. */ PetscFunctionBegin; PetscCall(PetscMalloc1(am + 1, &ci)); for (i = 0; i < am; i++) { const PetscInt anzi = ai[i + 1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */ const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */ a_rownnz = 0; for (k = 0; k < anzi; ++k) { a_rownnz += bi[acol[k] + 1] - bi[acol[k]]; if (a_rownnz > bn) { a_rownnz = bn; break; } } a_maxrownnz = PetscMax(a_maxrownnz, a_rownnz); } /* temporary work areas for merging rows */ PetscCall(PetscMalloc1(a_maxrownnz * 8, &workj_L1)); PetscCall(PetscMalloc1(a_maxrownnz * 8, &workj_L2)); PetscCall(PetscMalloc1(a_maxrownnz, &workj_L3)); /* This should be enough for almost all matrices. If not, memory is reallocated later. */ c_maxmem = 8 * (ai[am] + bi[bm]); /* Step 2: Populate pattern for C */ PetscCall(PetscMalloc1(c_maxmem, &cj)); ci_nnz = 0; ci[0] = 0; worki_L1[0] = 0; worki_L2[0] = 0; for (i = 0; i < am; i++) { const PetscInt anzi = ai[i + 1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */ const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */ rowsleft = anzi; inputcol_L1 = acol; L2_nnz = 0; L2_nrows = 1; /* Number of rows to be merged on Level 3. output of L3 already exists -> initial value 1 */ worki_L2[1] = 0; outputi_nnz = 0; /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem -> allocate more memory */ while (ci_nnz + a_maxrownnz > c_maxmem) { c_maxmem *= 2; ndouble++; PetscCall(PetscRealloc(sizeof(PetscInt) * c_maxmem, &cj)); } while (rowsleft) { L1_rowsleft = PetscMin(64, rowsleft); /* In the inner loop max 64 rows of B can be merged */ L1_nrows = 0; L1_nnz = 0; inputcol = inputcol_L1; inputi = bi; inputj = bj; /* The following macro is used to specialize for small rows in A. This helps with compiler unrolling, improving performance substantially. Input: inputj inputi inputcol bn Output: outputj outputi_nnz */ #define MatMatMultSymbolic_RowMergeMacro(ANNZ) \ do { \ window_min = bn; \ outputi_nnz = 0; \ for (k = 0; k < ANNZ; ++k) { \ brow_ptr[k] = inputj + inputi[inputcol[k]]; \ brow_end[k] = inputj + inputi[inputcol[k] + 1]; \ window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \ window_min = PetscMin(window[k], window_min); \ } \ while (window_min < bn) { \ outputj[outputi_nnz++] = window_min; \ /* advance front and compute new minimum */ \ old_window_min = window_min; \ window_min = bn; \ for (k = 0; k < ANNZ; ++k) { \ if (window[k] == old_window_min) { \ brow_ptr[k]++; \ window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \ } \ window_min = PetscMin(window[k], window_min); \ } \ } \ } while (0) /************** L E V E L 1 ***************/ /* Merge up to 8 rows of B to L1 work array*/ while (L1_rowsleft) { outputi_nnz = 0; if (anzi > 8) outputj = workj_L1 + L1_nnz; /* Level 1 rowmerge*/ else outputj = cj + ci_nnz; /* Merge directly to C */ switch (L1_rowsleft) { case 1: brow_ptr[0] = inputj + inputi[inputcol[0]]; brow_end[0] = inputj + inputi[inputcol[0] + 1]; for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */ inputcol += L1_rowsleft; rowsleft -= L1_rowsleft; L1_rowsleft = 0; break; case 2: MatMatMultSymbolic_RowMergeMacro(2); inputcol += L1_rowsleft; rowsleft -= L1_rowsleft; L1_rowsleft = 0; break; case 3: MatMatMultSymbolic_RowMergeMacro(3); inputcol += L1_rowsleft; rowsleft -= L1_rowsleft; L1_rowsleft = 0; break; case 4: MatMatMultSymbolic_RowMergeMacro(4); inputcol += L1_rowsleft; rowsleft -= L1_rowsleft; L1_rowsleft = 0; break; case 5: MatMatMultSymbolic_RowMergeMacro(5); inputcol += L1_rowsleft; rowsleft -= L1_rowsleft; L1_rowsleft = 0; break; case 6: MatMatMultSymbolic_RowMergeMacro(6); inputcol += L1_rowsleft; rowsleft -= L1_rowsleft; L1_rowsleft = 0; break; case 7: MatMatMultSymbolic_RowMergeMacro(7); inputcol += L1_rowsleft; rowsleft -= L1_rowsleft; L1_rowsleft = 0; break; default: MatMatMultSymbolic_RowMergeMacro(8); inputcol += 8; rowsleft -= 8; L1_rowsleft -= 8; break; } inputcol_L1 = inputcol; L1_nnz += outputi_nnz; worki_L1[++L1_nrows] = L1_nnz; } /********************** L E V E L 2 ************************/ /* Merge from L1 work array to either C or to L2 work array */ if (anzi > 8) { inputi = worki_L1; inputj = workj_L1; inputcol = workcol; outputi_nnz = 0; if (anzi <= 64) outputj = cj + ci_nnz; /* Merge from L1 work array to C */ else outputj = workj_L2 + L2_nnz; /* Merge from L1 work array to L2 work array */ switch (L1_nrows) { case 1: brow_ptr[0] = inputj + inputi[inputcol[0]]; brow_end[0] = inputj + inputi[inputcol[0] + 1]; for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */ break; case 2: MatMatMultSymbolic_RowMergeMacro(2); break; case 3: MatMatMultSymbolic_RowMergeMacro(3); break; case 4: MatMatMultSymbolic_RowMergeMacro(4); break; case 5: MatMatMultSymbolic_RowMergeMacro(5); break; case 6: MatMatMultSymbolic_RowMergeMacro(6); break; case 7: MatMatMultSymbolic_RowMergeMacro(7); break; case 8: MatMatMultSymbolic_RowMergeMacro(8); break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatMatMult logic error: Not merging 1-8 rows from L1 work array!"); } L2_nnz += outputi_nnz; worki_L2[++L2_nrows] = L2_nnz; /************************ L E V E L 3 **********************/ /* Merge from L2 work array to either C or to L2 work array */ if (anzi > 64 && (L2_nrows == 8 || rowsleft == 0)) { inputi = worki_L2; inputj = workj_L2; inputcol = workcol; outputi_nnz = 0; if (rowsleft) outputj = workj_L3; else outputj = cj + ci_nnz; switch (L2_nrows) { case 1: brow_ptr[0] = inputj + inputi[inputcol[0]]; brow_end[0] = inputj + inputi[inputcol[0] + 1]; for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */ break; case 2: MatMatMultSymbolic_RowMergeMacro(2); break; case 3: MatMatMultSymbolic_RowMergeMacro(3); break; case 4: MatMatMultSymbolic_RowMergeMacro(4); break; case 5: MatMatMultSymbolic_RowMergeMacro(5); break; case 6: MatMatMultSymbolic_RowMergeMacro(6); break; case 7: MatMatMultSymbolic_RowMergeMacro(7); break; case 8: MatMatMultSymbolic_RowMergeMacro(8); break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatMatMult logic error: Not merging 1-8 rows from L2 work array!"); } L2_nrows = 1; L2_nnz = outputi_nnz; worki_L2[1] = outputi_nnz; /* Copy to workj_L2 */ if (rowsleft) { for (k = 0; k < outputi_nnz; ++k) workj_L2[k] = outputj[k]; } } } } /* while (rowsleft) */ #undef MatMatMultSymbolic_RowMergeMacro /* terminate current row */ ci_nnz += outputi_nnz; ci[i + 1] = ci_nnz; } /* Step 3: Create the new symbolic matrix */ PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C)); PetscCall(MatSetBlockSizesFromMats(C, A, B)); /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ c = (Mat_SeqAIJ *)C->data; c->free_a = PETSC_TRUE; c->free_ij = PETSC_TRUE; c->nonew = 0; C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted; /* set MatInfo */ afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5; if (afill < 1.0) afill = 1.0; C->info.mallocs = ndouble; C->info.fill_ratio_given = fill; C->info.fill_ratio_needed = afill; #if defined(PETSC_USE_INFO) if (ci[am]) { PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill)); PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill)); } else { PetscCall(PetscInfo(C, "Empty matrix product\n")); } #endif /* Step 4: Free temporary work areas */ PetscCall(PetscFree(workj_L1)); PetscCall(PetscFree(workj_L2)); PetscCall(PetscFree(workj_L3)); PetscFunctionReturn(PETSC_SUCCESS); } /* concatenate unique entries and then sort */ PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(Mat A, Mat B, PetscReal fill, Mat C) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c; const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j; PetscInt *ci, *cj, bcol; PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N; PetscReal afill; PetscInt i, j, ndouble = 0; PetscSegBuffer seg, segrow; char *seen; PetscFunctionBegin; PetscCall(PetscMalloc1(am + 1, &ci)); ci[0] = 0; /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */ PetscCall(PetscSegBufferCreate(sizeof(PetscInt), (PetscInt)(fill * (ai[am] + bi[bm])), &seg)); PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 100, &segrow)); PetscCall(PetscCalloc1(bn, &seen)); /* Determine ci and cj */ for (i = 0; i < am; i++) { const PetscInt anzi = ai[i + 1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */ const PetscInt *acol = PetscSafePointerPlusOffset(aj, ai[i]); /* column indices of nonzero entries in this row */ PetscInt packlen = 0, *PETSC_RESTRICT crow; /* Pack segrow */ for (j = 0; j < anzi; j++) { PetscInt brow = acol[j], bjstart = bi[brow], bjend = bi[brow + 1], k; for (k = bjstart; k < bjend; k++) { bcol = bj[k]; if (!seen[bcol]) { /* new entry */ PetscInt *PETSC_RESTRICT slot; PetscCall(PetscSegBufferGetInts(segrow, 1, &slot)); *slot = bcol; seen[bcol] = 1; packlen++; } } } /* Check i-th diagonal entry */ if (C->force_diagonals && !seen[i]) { PetscInt *PETSC_RESTRICT slot; PetscCall(PetscSegBufferGetInts(segrow, 1, &slot)); *slot = i; seen[i] = 1; packlen++; } PetscCall(PetscSegBufferGetInts(seg, packlen, &crow)); PetscCall(PetscSegBufferExtractTo(segrow, crow)); PetscCall(PetscSortInt(packlen, crow)); ci[i + 1] = ci[i] + packlen; for (j = 0; j < packlen; j++) seen[crow[j]] = 0; } PetscCall(PetscSegBufferDestroy(&segrow)); PetscCall(PetscFree(seen)); /* Column indices are in the segmented buffer */ PetscCall(PetscSegBufferExtractAlloc(seg, &cj)); PetscCall(PetscSegBufferDestroy(&seg)); /* put together the new symbolic matrix */ PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C)); PetscCall(MatSetBlockSizesFromMats(C, A, B)); /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */ /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */ c = (Mat_SeqAIJ *)C->data; c->free_a = PETSC_TRUE; c->free_ij = PETSC_TRUE; c->nonew = 0; C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted; /* set MatInfo */ afill = (PetscReal)ci[am] / PetscMax(ai[am] + bi[bm], 1) + 1.e-5; if (afill < 1.0) afill = 1.0; C->info.mallocs = ndouble; C->info.fill_ratio_given = fill; C->info.fill_ratio_needed = afill; #if defined(PETSC_USE_INFO) if (ci[am]) { PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill)); PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill)); } else { PetscCall(PetscInfo(C, "Empty matrix product\n")); } #endif PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductCtxDestroy_SeqAIJ_MatMatMultTrans(PetscCtxRt data) { MatProductCtx_MatMatTransMult *abt = *(MatProductCtx_MatMatTransMult **)data; PetscFunctionBegin; PetscCall(MatTransposeColoringDestroy(&abt->matcoloring)); PetscCall(MatDestroy(&abt->Bt_den)); PetscCall(MatDestroy(&abt->ABt_den)); PetscCall(PetscFree(abt)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C) { Mat Bt; MatProductCtx_MatMatTransMult *abt; Mat_Product *product = C->product; char *alg; PetscFunctionBegin; PetscCheck(product, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct"); PetscCheck(!product->data, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Extra product struct not empty"); /* create symbolic Bt */ PetscCall(MatTransposeSymbolic(B, &Bt)); PetscCall(MatSetBlockSizes(Bt, A->cmap->bs, B->cmap->bs)); PetscCall(MatSetType(Bt, ((PetscObject)A)->type_name)); /* get symbolic C=A*Bt */ PetscCall(PetscStrallocpy(product->alg, &alg)); PetscCall(MatProductSetAlgorithm(C, "sorted")); /* set algorithm for C = A*Bt */ PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(A, Bt, fill, C)); PetscCall(MatProductSetAlgorithm(C, alg)); /* resume original algorithm for ABt product */ PetscCall(PetscFree(alg)); /* create a supporting struct for reuse intermediate dense matrices with matcoloring */ PetscCall(PetscNew(&abt)); product->data = abt; product->destroy = MatProductCtxDestroy_SeqAIJ_MatMatMultTrans; C->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ; abt->usecoloring = PETSC_FALSE; PetscCall(PetscStrcmp(product->alg, "color", &abt->usecoloring)); if (abt->usecoloring) { /* Create MatTransposeColoring from symbolic C=A*B^T */ MatTransposeColoring matcoloring; MatColoring coloring; ISColoring iscoloring; Mat Bt_dense, C_dense; /* inode causes memory problem */ PetscCall(MatSetOption(C, MAT_USE_INODES, PETSC_FALSE)); PetscCall(MatColoringCreate(C, &coloring)); PetscCall(MatColoringSetDistance(coloring, 2)); PetscCall(MatColoringSetType(coloring, MATCOLORINGSL)); PetscCall(MatColoringSetFromOptions(coloring)); PetscCall(MatColoringApply(coloring, &iscoloring)); PetscCall(MatColoringDestroy(&coloring)); PetscCall(MatTransposeColoringCreate(C, iscoloring, &matcoloring)); abt->matcoloring = matcoloring; PetscCall(ISColoringDestroy(&iscoloring)); /* Create Bt_dense and C_dense = A*Bt_dense */ PetscCall(MatCreate(PETSC_COMM_SELF, &Bt_dense)); PetscCall(MatSetSizes(Bt_dense, A->cmap->n, matcoloring->ncolors, A->cmap->n, matcoloring->ncolors)); PetscCall(MatSetType(Bt_dense, MATSEQDENSE)); PetscCall(MatSeqDenseSetPreallocation(Bt_dense, NULL)); Bt_dense->assembled = PETSC_TRUE; abt->Bt_den = Bt_dense; PetscCall(MatCreate(PETSC_COMM_SELF, &C_dense)); PetscCall(MatSetSizes(C_dense, A->rmap->n, matcoloring->ncolors, A->rmap->n, matcoloring->ncolors)); PetscCall(MatSetType(C_dense, MATSEQDENSE)); PetscCall(MatSeqDenseSetPreallocation(C_dense, NULL)); Bt_dense->assembled = PETSC_TRUE; abt->ABt_den = C_dense; #if defined(PETSC_USE_INFO) { Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; PetscCall(PetscInfo(C, "Use coloring of C=A*B^T; B^T: %" PetscInt_FMT " %" PetscInt_FMT ", Bt_dense: %" PetscInt_FMT ",%" PetscInt_FMT "; Cnz %" PetscInt_FMT " / (cm*ncolors %" PetscInt_FMT ") = %g\n", B->cmap->n, B->rmap->n, Bt_dense->rmap->n, Bt_dense->cmap->n, c->nz, A->rmap->n * matcoloring->ncolors, (double)(((PetscReal)c->nz) / ((PetscReal)(A->rmap->n * matcoloring->ncolors))))); } #endif } /* clean up */ PetscCall(MatDestroy(&Bt)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c = (Mat_SeqAIJ *)C->data; PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, anzi, bnzj, nexta, nextb, *acol, *bcol, brow; PetscInt cm = C->rmap->n, *ci = c->i, *cj = c->j, i, j, cnzi, *ccol; PetscLogDouble flops = 0.0; MatScalar *aa = a->a, *aval, *ba = b->a, *bval, *ca, *cval; MatProductCtx_MatMatTransMult *abt; Mat_Product *product = C->product; PetscFunctionBegin; PetscCheck(product, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct"); abt = (MatProductCtx_MatMatTransMult *)product->data; PetscCheck(abt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct"); /* clear old values in C */ if (!c->a) { PetscCall(PetscCalloc1(ci[cm] + 1, &ca)); c->a = ca; c->free_a = PETSC_TRUE; } else { ca = c->a; PetscCall(PetscArrayzero(ca, ci[cm] + 1)); } if (abt->usecoloring) { MatTransposeColoring matcoloring = abt->matcoloring; Mat Bt_dense, C_dense = abt->ABt_den; /* Get Bt_dense by Apply MatTransposeColoring to B */ Bt_dense = abt->Bt_den; PetscCall(MatTransColoringApplySpToDen(matcoloring, B, Bt_dense)); /* C_dense = A*Bt_dense */ PetscCall(MatMatMultNumeric_SeqAIJ_SeqDense(A, Bt_dense, C_dense)); /* Recover C from C_dense */ PetscCall(MatTransColoringApplyDenToSp(matcoloring, C_dense, C)); PetscFunctionReturn(PETSC_SUCCESS); } for (i = 0; i < cm; i++) { anzi = ai[i + 1] - ai[i]; acol = PetscSafePointerPlusOffset(aj, ai[i]); aval = PetscSafePointerPlusOffset(aa, ai[i]); cnzi = ci[i + 1] - ci[i]; ccol = PetscSafePointerPlusOffset(cj, ci[i]); cval = ca + ci[i]; for (j = 0; j < cnzi; j++) { brow = ccol[j]; bnzj = bi[brow + 1] - bi[brow]; bcol = bj + bi[brow]; bval = ba + bi[brow]; /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */ nexta = 0; nextb = 0; while (nexta < anzi && nextb < bnzj) { while (nexta < anzi && acol[nexta] < bcol[nextb]) nexta++; if (nexta == anzi) break; while (nextb < bnzj && acol[nexta] > bcol[nextb]) nextb++; if (nextb == bnzj) break; if (acol[nexta] == bcol[nextb]) { cval[j] += aval[nexta] * bval[nextb]; nexta++; nextb++; flops += 2; } } } } PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); PetscCall(PetscLogFlops(flops)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatProductCtxDestroy_SeqAIJ_MatTransMatMult(PetscCtxRt data) { MatProductCtx_MatTransMatMult *atb = *(MatProductCtx_MatTransMatMult **)data; PetscFunctionBegin; PetscCall(MatDestroy(&atb->At)); if (atb->destroy) PetscCall((*atb->destroy)(&atb->data)); PetscCall(PetscFree(atb)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C) { Mat At = NULL; Mat_Product *product = C->product; PetscBool flg, def, square; PetscFunctionBegin; MatCheckProduct(C, 4); square = (PetscBool)(A == B && A->symmetric == PETSC_BOOL3_TRUE); /* outerproduct */ PetscCall(PetscStrcmp(product->alg, "outerproduct", &flg)); if (flg) { /* create symbolic At */ if (!square) { PetscCall(MatTransposeSymbolic(A, &At)); PetscCall(MatSetBlockSizes(At, A->cmap->bs, B->cmap->bs)); PetscCall(MatSetType(At, ((PetscObject)A)->type_name)); } /* get symbolic C=At*B */ PetscCall(MatProductSetAlgorithm(C, "sorted")); PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(square ? A : At, B, fill, C)); /* clean up */ if (!square) PetscCall(MatDestroy(&At)); C->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ; /* outerproduct */ PetscCall(MatProductSetAlgorithm(C, "outerproduct")); PetscFunctionReturn(PETSC_SUCCESS); } /* matmatmult */ PetscCall(PetscStrcmp(product->alg, "default", &def)); PetscCall(PetscStrcmp(product->alg, "at*b", &flg)); if (flg || def) { MatProductCtx_MatTransMatMult *atb; PetscCheck(!product->data, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Extra product struct not empty"); PetscCall(PetscNew(&atb)); if (!square) PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &At)); PetscCall(MatProductSetAlgorithm(C, "sorted")); PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(square ? A : At, B, fill, C)); PetscCall(MatProductSetAlgorithm(C, "at*b")); product->data = atb; product->destroy = MatProductCtxDestroy_SeqAIJ_MatTransMatMult; atb->At = At; C->ops->mattransposemultnumeric = NULL; /* see MatProductNumeric_AtB_SeqAIJ_SeqAIJ */ PetscFunctionReturn(PETSC_SUCCESS); } SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Mat Product Algorithm is not supported"); } PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c = (Mat_SeqAIJ *)C->data; PetscInt am = A->rmap->n, anzi, *ai = a->i, *aj = a->j, *bi = b->i, *bj, bnzi, nextb; PetscInt cm = C->rmap->n, *ci = c->i, *cj = c->j, crow, *cjj, i, j, k; PetscLogDouble flops = 0.0; MatScalar *aa = a->a, *ba, *ca, *caj; PetscFunctionBegin; if (!c->a) { PetscCall(PetscCalloc1(ci[cm] + 1, &ca)); c->a = ca; c->free_a = PETSC_TRUE; } else { ca = c->a; PetscCall(PetscArrayzero(ca, ci[cm])); } /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */ for (i = 0; i < am; i++) { bj = b->j + bi[i]; ba = b->a + bi[i]; bnzi = bi[i + 1] - bi[i]; anzi = ai[i + 1] - ai[i]; for (j = 0; j < anzi; j++) { nextb = 0; crow = *aj++; cjj = cj + ci[crow]; caj = ca + ci[crow]; /* perform sparse axpy operation. Note cjj includes bj. */ for (k = 0; nextb < bnzi; k++) { if (cjj[k] == *(bj + nextb)) { /* ccol == bcol */ caj[k] += (*aa) * (*(ba + nextb)); nextb++; } } flops += 2 * bnzi; aa++; } } /* Assemble the final matrix and clean up */ PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); PetscCall(PetscLogFlops(flops)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A, Mat B, PetscReal fill, Mat C) { PetscFunctionBegin; PetscCall(MatMatMultSymbolic_SeqDense_SeqDense(A, B, 0.0, C)); C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; PetscFunctionReturn(PETSC_SUCCESS); } PETSC_INTERN PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A, Mat B, Mat C, const PetscBool add) { Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; PetscScalar *c, r1, r2, r3, r4, *c1, *c2, *c3, *c4; const PetscScalar *aa, *b, *b1, *b2, *b3, *b4, *av; const PetscInt *aj; PetscInt cm = C->rmap->n, cn = B->cmap->n, bm, am = A->rmap->n; PetscInt clda; PetscInt am4, bm4, col, i, j, n; PetscFunctionBegin; if (!cm || !cn) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(MatSeqAIJGetArrayRead(A, &av)); if (add) { PetscCall(MatDenseGetArray(C, &c)); } else { PetscCall(MatDenseGetArrayWrite(C, &c)); } PetscCall(MatDenseGetArrayRead(B, &b)); PetscCall(MatDenseGetLDA(B, &bm)); PetscCall(MatDenseGetLDA(C, &clda)); am4 = 4 * clda; bm4 = 4 * bm; if (b) { b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm; } else b1 = b2 = b3 = b4 = NULL; c1 = c; c2 = c1 + clda; c3 = c2 + clda; c4 = c3 + clda; for (col = 0; col < (cn / 4) * 4; col += 4) { /* over columns of C */ for (i = 0; i < am; i++) { /* over rows of A in those columns */ r1 = r2 = r3 = r4 = 0.0; n = a->i[i + 1] - a->i[i]; aj = PetscSafePointerPlusOffset(a->j, a->i[i]); aa = PetscSafePointerPlusOffset(av, a->i[i]); for (j = 0; j < n; j++) { const PetscScalar aatmp = aa[j]; const PetscInt ajtmp = aj[j]; r1 += aatmp * b1[ajtmp]; r2 += aatmp * b2[ajtmp]; r3 += aatmp * b3[ajtmp]; r4 += aatmp * b4[ajtmp]; } if (add) { c1[i] += r1; c2[i] += r2; c3[i] += r3; c4[i] += r4; } else { c1[i] = r1; c2[i] = r2; c3[i] = r3; c4[i] = r4; } } if (b) { b1 += bm4; b2 += bm4; b3 += bm4; b4 += bm4; } c1 += am4; c2 += am4; c3 += am4; c4 += am4; } /* process remaining columns */ if (col != cn) { PetscInt rc = cn - col; if (rc == 1) { for (i = 0; i < am; i++) { r1 = 0.0; n = a->i[i + 1] - a->i[i]; aj = PetscSafePointerPlusOffset(a->j, a->i[i]); aa = PetscSafePointerPlusOffset(av, a->i[i]); for (j = 0; j < n; j++) r1 += aa[j] * b1[aj[j]]; if (add) c1[i] += r1; else c1[i] = r1; } } else if (rc == 2) { for (i = 0; i < am; i++) { r1 = r2 = 0.0; n = a->i[i + 1] - a->i[i]; aj = PetscSafePointerPlusOffset(a->j, a->i[i]); aa = PetscSafePointerPlusOffset(av, a->i[i]); for (j = 0; j < n; j++) { const PetscScalar aatmp = aa[j]; const PetscInt ajtmp = aj[j]; r1 += aatmp * b1[ajtmp]; r2 += aatmp * b2[ajtmp]; } if (add) { c1[i] += r1; c2[i] += r2; } else { c1[i] = r1; c2[i] = r2; } } } else { for (i = 0; i < am; i++) { r1 = r2 = r3 = 0.0; n = a->i[i + 1] - a->i[i]; aj = PetscSafePointerPlusOffset(a->j, a->i[i]); aa = PetscSafePointerPlusOffset(av, a->i[i]); for (j = 0; j < n; j++) { const PetscScalar aatmp = aa[j]; const PetscInt ajtmp = aj[j]; r1 += aatmp * b1[ajtmp]; r2 += aatmp * b2[ajtmp]; r3 += aatmp * b3[ajtmp]; } if (add) { c1[i] += r1; c2[i] += r2; c3[i] += r3; } else { c1[i] = r1; c2[i] = r2; c3[i] = r3; } } } } PetscCall(PetscLogFlops(cn * (2.0 * a->nz))); if (add) { PetscCall(MatDenseRestoreArray(C, &c)); } else { PetscCall(MatDenseRestoreArrayWrite(C, &c)); } PetscCall(MatDenseRestoreArrayRead(B, &b)); PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A, Mat B, Mat C) { PetscFunctionBegin; PetscCheck(B->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number columns in A %" PetscInt_FMT " not equal rows in B %" PetscInt_FMT, A->cmap->n, B->rmap->n); PetscCheck(A->rmap->n == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number rows in C %" PetscInt_FMT " not equal rows in A %" PetscInt_FMT, C->rmap->n, A->rmap->n); PetscCheck(B->cmap->n == C->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number columns in B %" PetscInt_FMT " not equal columns in C %" PetscInt_FMT, B->cmap->n, C->cmap->n); PetscCall(MatMatMultNumericAdd_SeqAIJ_SeqDense(A, B, C, PETSC_FALSE)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_AB(Mat C) { PetscFunctionBegin; C->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJ_SeqDense; C->ops->productsymbolic = MatProductSymbolic_AB; PetscFunctionReturn(PETSC_SUCCESS); } PETSC_INTERN PetscErrorCode MatTMatTMultSymbolic_SeqAIJ_SeqDense(Mat, Mat, PetscReal, Mat); static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_AtB(Mat C) { PetscFunctionBegin; C->ops->transposematmultsymbolic = MatTMatTMultSymbolic_SeqAIJ_SeqDense; C->ops->productsymbolic = MatProductSymbolic_AtB; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_ABt(Mat C) { PetscFunctionBegin; C->ops->mattransposemultsymbolic = MatTMatTMultSymbolic_SeqAIJ_SeqDense; C->ops->productsymbolic = MatProductSymbolic_ABt; PetscFunctionReturn(PETSC_SUCCESS); } PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense(Mat C) { Mat_Product *product = C->product; PetscFunctionBegin; switch (product->type) { case MATPRODUCT_AB: PetscCall(MatProductSetFromOptions_SeqAIJ_SeqDense_AB(C)); break; case MATPRODUCT_AtB: PetscCall(MatProductSetFromOptions_SeqAIJ_SeqDense_AtB(C)); break; case MATPRODUCT_ABt: PetscCall(MatProductSetFromOptions_SeqAIJ_SeqDense_ABt(C)); break; default: break; } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSetFromOptions_SeqXBAIJ_SeqDense_AB(Mat C) { Mat_Product *product = C->product; Mat A = product->A; PetscBool baij; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQBAIJ, &baij)); if (!baij) { /* A is seqsbaij */ PetscBool sbaij; PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &sbaij)); PetscCheck(sbaij, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Mat must be either seqbaij or seqsbaij format"); C->ops->matmultsymbolic = MatMatMultSymbolic_SeqSBAIJ_SeqDense; } else { /* A is seqbaij */ C->ops->matmultsymbolic = MatMatMultSymbolic_SeqBAIJ_SeqDense; } C->ops->productsymbolic = MatProductSymbolic_AB; PetscFunctionReturn(PETSC_SUCCESS); } PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqXBAIJ_SeqDense(Mat C) { Mat_Product *product = C->product; PetscFunctionBegin; MatCheckProduct(C, 1); PetscCheck(product->A, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing A"); if (product->type == MATPRODUCT_AB || (product->type == MATPRODUCT_AtB && product->A->symmetric == PETSC_BOOL3_TRUE)) PetscCall(MatProductSetFromOptions_SeqXBAIJ_SeqDense_AB(C)); else if (product->type == MATPRODUCT_AtB) { PetscBool flg; PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATSEQBAIJ, &flg)); if (flg) { C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqBAIJ_SeqDense; C->ops->productsymbolic = MatProductSymbolic_AtB; } } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSetFromOptions_SeqDense_SeqAIJ_AB(Mat C) { PetscFunctionBegin; C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqAIJ; C->ops->productsymbolic = MatProductSymbolic_AB; PetscFunctionReturn(PETSC_SUCCESS); } PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense_SeqAIJ(Mat C) { Mat_Product *product = C->product; PetscFunctionBegin; if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_SeqDense_SeqAIJ_AB(C)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatTransColoringApplySpToDen_SeqAIJ(MatTransposeColoring coloring, Mat B, Mat Btdense) { Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; Mat_SeqDense *btdense = (Mat_SeqDense *)Btdense->data; PetscInt *bi = b->i, *bj = b->j; PetscInt m = Btdense->rmap->n, n = Btdense->cmap->n, j, k, l, col, anz, *btcol, brow, ncolumns; MatScalar *btval, *btval_den, *ba = b->a; PetscInt *columns = coloring->columns, *colorforcol = coloring->colorforcol, ncolors = coloring->ncolors; PetscFunctionBegin; btval_den = btdense->v; PetscCall(PetscArrayzero(btval_den, m * n)); for (k = 0; k < ncolors; k++) { ncolumns = coloring->ncolumns[k]; for (l = 0; l < ncolumns; l++) { /* insert a row of B to a column of Btdense */ col = *(columns + colorforcol[k] + l); btcol = bj + bi[col]; btval = ba + bi[col]; anz = bi[col + 1] - bi[col]; for (j = 0; j < anz; j++) { brow = btcol[j]; btval_den[brow] = btval[j]; } } btval_den += m; } PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring, Mat Cden, Mat Csp) { Mat_SeqAIJ *csp = (Mat_SeqAIJ *)Csp->data; const PetscScalar *ca_den, *ca_den_ptr; PetscScalar *ca = csp->a; PetscInt k, l, m = Cden->rmap->n, ncolors = matcoloring->ncolors; PetscInt brows = matcoloring->brows, *den2sp = matcoloring->den2sp; PetscInt nrows, *row, *idx; PetscInt *rows = matcoloring->rows, *colorforrow = matcoloring->colorforrow; PetscFunctionBegin; PetscCall(MatDenseGetArrayRead(Cden, &ca_den)); if (brows > 0) { PetscInt *lstart, row_end, row_start; lstart = matcoloring->lstart; PetscCall(PetscArrayzero(lstart, ncolors)); row_end = brows; if (row_end > m) row_end = m; for (row_start = 0; row_start < m; row_start += brows) { /* loop over row blocks of Csp */ ca_den_ptr = ca_den; for (k = 0; k < ncolors; k++) { /* loop over colors (columns of Cden) */ nrows = matcoloring->nrows[k]; row = rows + colorforrow[k]; idx = den2sp + colorforrow[k]; for (l = lstart[k]; l < nrows; l++) { if (row[l] >= row_end) { lstart[k] = l; break; } else { ca[idx[l]] = ca_den_ptr[row[l]]; } } ca_den_ptr += m; } row_end += brows; if (row_end > m) row_end = m; } } else { /* non-blocked impl: loop over columns of Csp - slow if Csp is large */ ca_den_ptr = ca_den; for (k = 0; k < ncolors; k++) { nrows = matcoloring->nrows[k]; row = rows + colorforrow[k]; idx = den2sp + colorforrow[k]; for (l = 0; l < nrows; l++) ca[idx[l]] = ca_den_ptr[row[l]]; ca_den_ptr += m; } } PetscCall(MatDenseRestoreArrayRead(Cden, &ca_den)); #if defined(PETSC_USE_INFO) if (matcoloring->brows > 0) { PetscCall(PetscInfo(Csp, "Loop over %" PetscInt_FMT " row blocks for den2sp\n", brows)); } else { PetscCall(PetscInfo(Csp, "Loop over colors/columns of Cden, inefficient for large sparse matrix product \n")); } #endif PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat, ISColoring iscoloring, MatTransposeColoring c) { PetscInt i, n, nrows, Nbs, j, k, m, ncols, col, cm; const PetscInt *is, *ci, *cj, *row_idx; PetscInt nis = iscoloring->n, *rowhit, bs = 1; IS *isa; Mat_SeqAIJ *csp = (Mat_SeqAIJ *)mat->data; PetscInt *colorforrow, *rows, *rows_i, *idxhit, *spidx, *den2sp, *den2sp_i; PetscInt *colorforcol, *columns, *columns_i, brows; PetscBool flg; PetscFunctionBegin; PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, PETSC_IGNORE, &isa)); /* bs > 1 is not being tested yet! */ Nbs = mat->cmap->N / bs; c->M = mat->rmap->N / bs; /* set total rows, columns and local rows */ c->N = Nbs; c->m = c->M; c->rstart = 0; c->brows = 100; c->ncolors = nis; PetscCall(PetscMalloc3(nis, &c->ncolumns, nis, &c->nrows, nis + 1, &colorforrow)); PetscCall(PetscMalloc1(csp->nz + 1, &rows)); PetscCall(PetscMalloc1(csp->nz + 1, &den2sp)); brows = c->brows; PetscCall(PetscOptionsGetInt(NULL, NULL, "-matden2sp_brows", &brows, &flg)); if (flg) c->brows = brows; if (brows > 0) PetscCall(PetscMalloc1(nis + 1, &c->lstart)); colorforrow[0] = 0; rows_i = rows; den2sp_i = den2sp; PetscCall(PetscMalloc1(nis + 1, &colorforcol)); PetscCall(PetscMalloc1(Nbs + 1, &columns)); colorforcol[0] = 0; columns_i = columns; /* get column-wise storage of mat */ PetscCall(MatGetColumnIJ_SeqAIJ_Color(mat, 0, PETSC_FALSE, PETSC_FALSE, &ncols, &ci, &cj, &spidx, NULL)); cm = c->m; PetscCall(PetscMalloc1(cm + 1, &rowhit)); PetscCall(PetscMalloc1(cm + 1, &idxhit)); for (i = 0; i < nis; i++) { /* loop over color */ PetscCall(ISGetLocalSize(isa[i], &n)); PetscCall(ISGetIndices(isa[i], &is)); c->ncolumns[i] = n; if (n) PetscCall(PetscArraycpy(columns_i, is, n)); colorforcol[i + 1] = colorforcol[i] + n; columns_i += n; /* fast, crude version requires O(N*N) work */ PetscCall(PetscArrayzero(rowhit, cm)); for (j = 0; j < n; j++) { /* loop over columns*/ col = is[j]; row_idx = cj + ci[col]; m = ci[col + 1] - ci[col]; for (k = 0; k < m; k++) { /* loop over columns marking them in rowhit */ idxhit[*row_idx] = spidx[ci[col] + k]; rowhit[*row_idx++] = col + 1; } } /* count the number of hits */ nrows = 0; for (j = 0; j < cm; j++) { if (rowhit[j]) nrows++; } c->nrows[i] = nrows; colorforrow[i + 1] = colorforrow[i] + nrows; nrows = 0; for (j = 0; j < cm; j++) { /* loop over rows */ if (rowhit[j]) { rows_i[nrows] = j; den2sp_i[nrows] = idxhit[j]; nrows++; } } den2sp_i += nrows; PetscCall(ISRestoreIndices(isa[i], &is)); rows_i += nrows; } PetscCall(MatRestoreColumnIJ_SeqAIJ_Color(mat, 0, PETSC_FALSE, PETSC_FALSE, &ncols, &ci, &cj, &spidx, NULL)); PetscCall(PetscFree(rowhit)); PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &isa)); PetscCheck(csp->nz == colorforrow[nis], PETSC_COMM_SELF, PETSC_ERR_PLIB, "csp->nz %" PetscInt_FMT " != colorforrow[nis] %" PetscInt_FMT, csp->nz, colorforrow[nis]); c->colorforrow = colorforrow; c->rows = rows; c->den2sp = den2sp; c->colorforcol = colorforcol; c->columns = columns; PetscCall(PetscFree(idxhit)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductNumeric_AtB_SeqAIJ_SeqAIJ(Mat C) { Mat_Product *product = C->product; Mat A = product->A, B = product->B; PetscFunctionBegin; if (C->ops->mattransposemultnumeric) { /* Alg: "outerproduct" */ PetscCall((*C->ops->mattransposemultnumeric)(A, B, C)); } else { /* Alg: "matmatmult" -- C = At*B */ MatProductCtx_MatTransMatMult *atb = (MatProductCtx_MatTransMatMult *)product->data; PetscCheck(atb, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct"); if (atb->At) { /* At is computed in MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(); user may have called MatProductReplaceMats() to get this A=product->A */ PetscCall(MatTransposeSetPrecursor(A, atb->At)); PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &atb->At)); } PetscCall(MatMatMultNumeric_SeqAIJ_SeqAIJ(atb->At ? atb->At : A, B, C)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSymbolic_AtB_SeqAIJ_SeqAIJ(Mat C) { Mat_Product *product = C->product; Mat A = product->A, B = product->B; PetscReal fill = product->fill; PetscFunctionBegin; PetscCall(MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A, B, fill, C)); C->ops->productnumeric = MatProductNumeric_AtB_SeqAIJ_SeqAIJ; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSetFromOptions_SeqAIJ_AB(Mat C) { Mat_Product *product = C->product; PetscInt alg = 0; /* default algorithm */ PetscBool flg = PETSC_FALSE; #if !defined(PETSC_HAVE_HYPRE) const char *algTypes[7] = {"sorted", "scalable", "scalable_fast", "heap", "btheap", "llcondensed", "rowmerge"}; PetscInt nalg = 7; #else const char *algTypes[8] = {"sorted", "scalable", "scalable_fast", "heap", "btheap", "llcondensed", "rowmerge", "hypre"}; PetscInt nalg = 8; #endif PetscFunctionBegin; /* Set default algorithm */ PetscCall(PetscStrcmp(C->product->alg, "default", &flg)); if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); /* Get runtime option */ if (product->api_user) { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatMult", "Mat"); PetscCall(PetscOptionsEList("-matmatmult_via", "Algorithmic approach", "MatMatMult", algTypes, nalg, algTypes[0], &alg, &flg)); PetscOptionsEnd(); } else { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_AB", "Mat"); PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_AB", algTypes, nalg, algTypes[0], &alg, &flg)); PetscOptionsEnd(); } if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); C->ops->productsymbolic = MatProductSymbolic_AB; C->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJ_SeqAIJ; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSetFromOptions_SeqAIJ_AtB(Mat C) { Mat_Product *product = C->product; PetscInt alg = 0; /* default algorithm */ PetscBool flg = PETSC_FALSE; const char *algTypes[3] = {"default", "at*b", "outerproduct"}; PetscInt nalg = 3; PetscFunctionBegin; /* Get runtime option */ if (product->api_user) { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatTransposeMatMult", "Mat"); PetscCall(PetscOptionsEList("-mattransposematmult_via", "Algorithmic approach", "MatTransposeMatMult", algTypes, nalg, algTypes[alg], &alg, &flg)); PetscOptionsEnd(); } else { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_AtB", "Mat"); PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_AtB", algTypes, nalg, algTypes[alg], &alg, &flg)); PetscOptionsEnd(); } if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); C->ops->productsymbolic = MatProductSymbolic_AtB_SeqAIJ_SeqAIJ; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSetFromOptions_SeqAIJ_ABt(Mat C) { Mat_Product *product = C->product; PetscInt alg = 0; /* default algorithm */ PetscBool flg = PETSC_FALSE; const char *algTypes[2] = {"default", "color"}; PetscInt nalg = 2; PetscFunctionBegin; /* Set default algorithm */ PetscCall(PetscStrcmp(C->product->alg, "default", &flg)); if (!flg) { alg = 1; PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); } /* Get runtime option */ if (product->api_user) { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatTransposeMult", "Mat"); PetscCall(PetscOptionsEList("-matmattransmult_via", "Algorithmic approach", "MatMatTransposeMult", algTypes, nalg, algTypes[alg], &alg, &flg)); PetscOptionsEnd(); } else { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_ABt", "Mat"); PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_ABt", algTypes, nalg, algTypes[alg], &alg, &flg)); PetscOptionsEnd(); } if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ; C->ops->productsymbolic = MatProductSymbolic_ABt; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSetFromOptions_SeqAIJ_PtAP(Mat C) { Mat_Product *product = C->product; PetscBool flg = PETSC_FALSE; PetscInt alg = 0; /* default algorithm -- alg=1 should be default!!! */ #if !defined(PETSC_HAVE_HYPRE) const char *algTypes[2] = {"scalable", "rap"}; PetscInt nalg = 2; #else const char *algTypes[3] = {"scalable", "rap", "hypre"}; PetscInt nalg = 3; #endif PetscFunctionBegin; /* Set default algorithm */ PetscCall(PetscStrcmp(product->alg, "default", &flg)); if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); /* Get runtime option */ if (product->api_user) { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP", "Mat"); PetscCall(PetscOptionsEList("-matptap_via", "Algorithmic approach", "MatPtAP", algTypes, nalg, algTypes[0], &alg, &flg)); PetscOptionsEnd(); } else { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP", "Mat"); PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_PtAP", algTypes, nalg, algTypes[0], &alg, &flg)); PetscOptionsEnd(); } if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); C->ops->productsymbolic = MatProductSymbolic_PtAP_SeqAIJ_SeqAIJ; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode MatProductSetFromOptions_SeqAIJ_RARt(Mat C) { Mat_Product *product = C->product; PetscBool flg = PETSC_FALSE; PetscInt alg = 0; /* default algorithm */ const char *algTypes[3] = {"r*a*rt", "r*art", "coloring_rart"}; PetscInt nalg = 3; PetscFunctionBegin; /* Set default algorithm */ PetscCall(PetscStrcmp(product->alg, "default", &flg)); if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); /* Get runtime option */ if (product->api_user) { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatRARt", "Mat"); PetscCall(PetscOptionsEList("-matrart_via", "Algorithmic approach", "MatRARt", algTypes, nalg, algTypes[0], &alg, &flg)); PetscOptionsEnd(); } else { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_RARt", "Mat"); PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_RARt", algTypes, nalg, algTypes[0], &alg, &flg)); PetscOptionsEnd(); } if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); C->ops->productsymbolic = MatProductSymbolic_RARt_SeqAIJ_SeqAIJ; PetscFunctionReturn(PETSC_SUCCESS); } /* ABC = A*B*C = A*(B*C); ABC's algorithm must be chosen from AB's algorithm */ static PetscErrorCode MatProductSetFromOptions_SeqAIJ_ABC(Mat C) { Mat_Product *product = C->product; PetscInt alg = 0; /* default algorithm */ PetscBool flg = PETSC_FALSE; const char *algTypes[7] = {"sorted", "scalable", "scalable_fast", "heap", "btheap", "llcondensed", "rowmerge"}; PetscInt nalg = 7; PetscFunctionBegin; /* Set default algorithm */ PetscCall(PetscStrcmp(product->alg, "default", &flg)); if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); /* Get runtime option */ if (product->api_user) { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatMatMult", "Mat"); PetscCall(PetscOptionsEList("-matmatmatmult_via", "Algorithmic approach", "MatMatMatMult", algTypes, nalg, algTypes[alg], &alg, &flg)); PetscOptionsEnd(); } else { PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_ABC", "Mat"); PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_ABC", algTypes, nalg, algTypes[alg], &alg, &flg)); PetscOptionsEnd(); } if (flg) PetscCall(MatProductSetAlgorithm(C, algTypes[alg])); C->ops->matmatmultsymbolic = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ; C->ops->productsymbolic = MatProductSymbolic_ABC; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode MatProductSetFromOptions_SeqAIJ(Mat C) { Mat_Product *product = C->product; PetscFunctionBegin; switch (product->type) { case MATPRODUCT_AB: PetscCall(MatProductSetFromOptions_SeqAIJ_AB(C)); break; case MATPRODUCT_AtB: PetscCall(MatProductSetFromOptions_SeqAIJ_AtB(C)); break; case MATPRODUCT_ABt: PetscCall(MatProductSetFromOptions_SeqAIJ_ABt(C)); break; case MATPRODUCT_PtAP: PetscCall(MatProductSetFromOptions_SeqAIJ_PtAP(C)); break; case MATPRODUCT_RARt: PetscCall(MatProductSetFromOptions_SeqAIJ_RARt(C)); break; case MATPRODUCT_ABC: PetscCall(MatProductSetFromOptions_SeqAIJ_ABC(C)); break; default: break; } PetscFunctionReturn(PETSC_SUCCESS); }