1 2 /* 3 Defines a matrix-vector product for the MATSEQAIJCRL matrix class. 4 This class is derived from the MATSEQAIJ class and retains the 5 compressed row storage (aka Yale sparse matrix format) but augments 6 it with a column oriented storage that is more efficient for 7 matrix vector products on Vector machines. 8 9 CRL stands for constant row length (that is the same number of columns 10 is kept (padded with zeros) for each row of the sparse matrix. 11 */ 12 #include <../src/mat/impls/aij/seq/crl/crl.h> 13 14 PetscErrorCode MatDestroy_SeqAIJCRL(Mat A) { 15 Mat_AIJCRL *aijcrl = (Mat_AIJCRL *)A->spptr; 16 17 PetscFunctionBegin; 18 /* Free everything in the Mat_AIJCRL data structure. */ 19 if (aijcrl) PetscCall(PetscFree2(aijcrl->acols, aijcrl->icols)); 20 PetscCall(PetscFree(A->spptr)); 21 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSEQAIJ)); 22 PetscCall(MatDestroy_SeqAIJ(A)); 23 PetscFunctionReturn(0); 24 } 25 26 PetscErrorCode MatDuplicate_AIJCRL(Mat A, MatDuplicateOption op, Mat *M) { 27 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot duplicate AIJCRL matrices yet"); 28 } 29 30 PetscErrorCode MatSeqAIJCRL_create_aijcrl(Mat A) { 31 Mat_SeqAIJ *a = (Mat_SeqAIJ *)(A)->data; 32 Mat_AIJCRL *aijcrl = (Mat_AIJCRL *)A->spptr; 33 PetscInt m = A->rmap->n; /* Number of rows in the matrix. */ 34 PetscInt *aj = a->j; /* From the CSR representation; points to the beginning of each row. */ 35 PetscInt i, j, rmax = a->rmax, *icols, *ilen = a->ilen; 36 MatScalar *aa = a->a; 37 PetscScalar *acols; 38 39 PetscFunctionBegin; 40 aijcrl->nz = a->nz; 41 aijcrl->m = A->rmap->n; 42 aijcrl->rmax = rmax; 43 44 PetscCall(PetscFree2(aijcrl->acols, aijcrl->icols)); 45 PetscCall(PetscMalloc2(rmax * m, &aijcrl->acols, rmax * m, &aijcrl->icols)); 46 acols = aijcrl->acols; 47 icols = aijcrl->icols; 48 for (i = 0; i < m; i++) { 49 for (j = 0; j < ilen[i]; j++) { 50 acols[j * m + i] = *aa++; 51 icols[j * m + i] = *aj++; 52 } 53 for (; j < rmax; j++) { /* empty column entries */ 54 acols[j * m + i] = 0.0; 55 icols[j * m + i] = (j) ? icols[(j - 1) * m + i] : 0; /* handle case where row is EMPTY */ 56 } 57 } 58 PetscCall(PetscInfo(A, "Percentage of 0's introduced for vectorized multiply %g. Rmax= %" PetscInt_FMT "\n", 1.0 - ((double)a->nz) / ((double)(rmax * m)), rmax)); 59 PetscFunctionReturn(0); 60 } 61 62 PetscErrorCode MatAssemblyEnd_SeqAIJCRL(Mat A, MatAssemblyType mode) { 63 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 64 65 PetscFunctionBegin; 66 a->inode.use = PETSC_FALSE; 67 68 PetscCall(MatAssemblyEnd_SeqAIJ(A, mode)); 69 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 70 71 /* Now calculate the permutation and grouping information. */ 72 PetscCall(MatSeqAIJCRL_create_aijcrl(A)); 73 PetscFunctionReturn(0); 74 } 75 76 #include <../src/mat/impls/aij/seq/crl/ftn-kernels/fmultcrl.h> 77 78 /* 79 Shared by both sequential and parallel versions of CRL matrix: MATMPIAIJCRL and MATSEQAIJCRL 80 - the scatter is used only in the parallel version 81 82 */ 83 PetscErrorCode MatMult_AIJCRL(Mat A, Vec xx, Vec yy) { 84 Mat_AIJCRL *aijcrl = (Mat_AIJCRL *)A->spptr; 85 PetscInt m = aijcrl->m; /* Number of rows in the matrix. */ 86 PetscInt rmax = aijcrl->rmax, *icols = aijcrl->icols; 87 PetscScalar *acols = aijcrl->acols; 88 PetscScalar *y; 89 const PetscScalar *x; 90 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTCRL) 91 PetscInt i, j, ii; 92 #endif 93 94 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 95 #pragma disjoint(*x, *y, *aa) 96 #endif 97 98 PetscFunctionBegin; 99 if (aijcrl->xscat) { 100 PetscCall(VecCopy(xx, aijcrl->xwork)); 101 /* get remote values needed for local part of multiply */ 102 PetscCall(VecScatterBegin(aijcrl->xscat, xx, aijcrl->fwork, INSERT_VALUES, SCATTER_FORWARD)); 103 PetscCall(VecScatterEnd(aijcrl->xscat, xx, aijcrl->fwork, INSERT_VALUES, SCATTER_FORWARD)); 104 xx = aijcrl->xwork; 105 } 106 107 PetscCall(VecGetArrayRead(xx, &x)); 108 PetscCall(VecGetArray(yy, &y)); 109 110 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTCRL) 111 fortranmultcrl_(&m, &rmax, x, y, icols, acols); 112 #else 113 114 /* first column */ 115 for (j = 0; j < m; j++) y[j] = acols[j] * x[icols[j]]; 116 117 /* other columns */ 118 #if defined(PETSC_HAVE_CRAY_VECTOR) 119 #pragma _CRI preferstream 120 #endif 121 for (i = 1; i < rmax; i++) { 122 ii = i * m; 123 #if defined(PETSC_HAVE_CRAY_VECTOR) 124 #pragma _CRI prefervector 125 #endif 126 for (j = 0; j < m; j++) y[j] = y[j] + acols[ii + j] * x[icols[ii + j]]; 127 } 128 #if defined(PETSC_HAVE_CRAY_VECTOR) 129 #pragma _CRI ivdep 130 #endif 131 132 #endif 133 PetscCall(PetscLogFlops(2.0 * aijcrl->nz - m)); 134 PetscCall(VecRestoreArrayRead(xx, &x)); 135 PetscCall(VecRestoreArray(yy, &y)); 136 PetscFunctionReturn(0); 137 } 138 139 /* MatConvert_SeqAIJ_SeqAIJCRL converts a SeqAIJ matrix into a 140 * SeqAIJCRL matrix. This routine is called by the MatCreate_SeqAIJCRL() 141 * routine, but can also be used to convert an assembled SeqAIJ matrix 142 * into a SeqAIJCRL one. */ 143 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat A, MatType type, MatReuse reuse, Mat *newmat) { 144 Mat B = *newmat; 145 Mat_AIJCRL *aijcrl; 146 PetscBool sametype; 147 148 PetscFunctionBegin; 149 if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B)); 150 PetscCall(PetscObjectTypeCompare((PetscObject)A, type, &sametype)); 151 if (sametype) PetscFunctionReturn(0); 152 153 PetscCall(PetscNewLog(B, &aijcrl)); 154 B->spptr = (void *)aijcrl; 155 156 /* Set function pointers for methods that we inherit from AIJ but override. */ 157 B->ops->duplicate = MatDuplicate_AIJCRL; 158 B->ops->assemblyend = MatAssemblyEnd_SeqAIJCRL; 159 B->ops->destroy = MatDestroy_SeqAIJCRL; 160 B->ops->mult = MatMult_AIJCRL; 161 162 /* If A has already been assembled, compute the permutation. */ 163 if (A->assembled) PetscCall(MatSeqAIJCRL_create_aijcrl(B)); 164 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJCRL)); 165 *newmat = B; 166 PetscFunctionReturn(0); 167 } 168 169 /*@C 170 MatCreateSeqAIJCRL - Creates a sparse matrix of type `MATSEQAIJCRL`. 171 This type inherits from `MATSEQAIJ`, but stores some additional 172 information that is used to allow better vectorization of 173 the matrix-vector product. At the cost of increased storage, the `MATSEQAIJ` formatted 174 matrix can be copied to a format in which pieces of the matrix are 175 stored in ELLPACK format, allowing the vectorized matrix multiply 176 routine to use stride-1 memory accesses. As with the `MATSEQAIJ` type, it is 177 important to preallocate matrix storage in order to get good assembly 178 performance. 179 180 Collective 181 182 Input Parameters: 183 + comm - MPI communicator, set to `PETSC_COMM_SELF` 184 . m - number of rows 185 . n - number of columns 186 . nz - number of nonzeros per row (same for all rows) 187 - nnz - array containing the number of nonzeros in the various rows 188 (possibly different for each row) or NULL 189 190 Output Parameter: 191 . A - the matrix 192 193 Note: 194 If nnz is given then nz is ignored 195 196 Level: intermediate 197 198 .seealso: `MatCreate()`, `MatCreateMPIAIJPERM()`, `MatSetValues()` 199 @*/ 200 PetscErrorCode MatCreateSeqAIJCRL(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A) { 201 PetscFunctionBegin; 202 PetscCall(MatCreate(comm, A)); 203 PetscCall(MatSetSizes(*A, m, n, m, n)); 204 PetscCall(MatSetType(*A, MATSEQAIJCRL)); 205 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A, nz, nnz)); 206 PetscFunctionReturn(0); 207 } 208 209 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJCRL(Mat A) { 210 PetscFunctionBegin; 211 PetscCall(MatSetType(A, MATSEQAIJ)); 212 PetscCall(MatConvert_SeqAIJ_SeqAIJCRL(A, MATSEQAIJCRL, MAT_INPLACE_MATRIX, &A)); 213 PetscFunctionReturn(0); 214 } 215