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