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