1 2 /* 3 Defines a matrix-vector product for the MATMPIAIJCRL matrix class. 4 This class is derived from the MATMPIAIJ 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 See src/mat/impls/aij/seq/crl/crl.c for the sequential version 13 */ 14 15 #include <../src/mat/impls/aij/mpi/mpiaij.h> 16 #include <../src/mat/impls/aij/seq/crl/crl.h> 17 18 PetscErrorCode MatDestroy_MPIAIJCRL(Mat A) 19 { 20 Mat_AIJCRL *aijcrl = (Mat_AIJCRL*) A->spptr; 21 22 PetscFunctionBegin; 23 if (aijcrl) { 24 PetscCall(PetscFree2(aijcrl->acols,aijcrl->icols)); 25 PetscCall(VecDestroy(&aijcrl->fwork)); 26 PetscCall(VecDestroy(&aijcrl->xwork)); 27 PetscCall(PetscFree(aijcrl->array)); 28 } 29 PetscCall(PetscFree(A->spptr)); 30 31 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATMPIAIJ)); 32 PetscCall(MatDestroy_MPIAIJ(A)); 33 PetscFunctionReturn(0); 34 } 35 36 PetscErrorCode MatMPIAIJCRL_create_aijcrl(Mat A) 37 { 38 Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A)->data; 39 Mat_SeqAIJ *Aij = (Mat_SeqAIJ*)(a->A->data), *Bij = (Mat_SeqAIJ*)(a->B->data); 40 Mat_AIJCRL *aijcrl = (Mat_AIJCRL*) A->spptr; 41 PetscInt m = A->rmap->n; /* Number of rows in the matrix. */ 42 PetscInt nd = a->A->cmap->n; /* number of columns in diagonal portion */ 43 PetscInt *aj = Aij->j,*bj = Bij->j; /* From the CSR representation; points to the beginning of each row. */ 44 PetscInt i, j,rmax = 0,*icols, *ailen = Aij->ilen, *bilen = Bij->ilen; 45 PetscScalar *aa = Aij->a,*ba = Bij->a,*acols,*array; 46 47 PetscFunctionBegin; 48 /* determine the row with the most columns */ 49 for (i=0; i<m; i++) { 50 rmax = PetscMax(rmax,ailen[i]+bilen[i]); 51 } 52 aijcrl->nz = Aij->nz+Bij->nz; 53 aijcrl->m = A->rmap->n; 54 aijcrl->rmax = rmax; 55 56 PetscCall(PetscFree2(aijcrl->acols,aijcrl->icols)); 57 PetscCall(PetscMalloc2(rmax*m,&aijcrl->acols,rmax*m,&aijcrl->icols)); 58 acols = aijcrl->acols; 59 icols = aijcrl->icols; 60 for (i=0; i<m; i++) { 61 for (j=0; j<ailen[i]; j++) { 62 acols[j*m+i] = *aa++; 63 icols[j*m+i] = *aj++; 64 } 65 for (; j<ailen[i]+bilen[i]; j++) { 66 acols[j*m+i] = *ba++; 67 icols[j*m+i] = nd + *bj++; 68 } 69 for (; j<rmax; j++) { /* empty column entries */ 70 acols[j*m+i] = 0.0; 71 icols[j*m+i] = (j) ? icols[(j-1)*m+i] : 0; /* handle case where row is EMPTY */ 72 } 73 } 74 PetscCall(PetscInfo(A,"Percentage of 0's introduced for vectorized multiply %g\n",1.0-((double)(aijcrl->nz))/((double)(rmax*m)))); 75 76 PetscCall(PetscFree(aijcrl->array)); 77 PetscCall(PetscMalloc1(a->B->cmap->n+nd,&array)); 78 /* xwork array is actually B->n+nd long, but we define xwork this length so can copy into it */ 79 PetscCall(VecDestroy(&aijcrl->xwork)); 80 PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,nd,PETSC_DECIDE,array,&aijcrl->xwork)); 81 PetscCall(VecDestroy(&aijcrl->fwork)); 82 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->B->cmap->n,array+nd,&aijcrl->fwork)); 83 84 aijcrl->array = array; 85 aijcrl->xscat = a->Mvctx; 86 PetscFunctionReturn(0); 87 } 88 89 PetscErrorCode MatAssemblyEnd_MPIAIJCRL(Mat A, MatAssemblyType mode) 90 { 91 Mat_MPIAIJ *a = (Mat_MPIAIJ*)A->data; 92 Mat_SeqAIJ *Aij = (Mat_SeqAIJ*)(a->A->data), *Bij = (Mat_SeqAIJ*)(a->A->data); 93 94 PetscFunctionBegin; 95 Aij->inode.use = PETSC_FALSE; 96 Bij->inode.use = PETSC_FALSE; 97 98 PetscCall(MatAssemblyEnd_MPIAIJ(A,mode)); 99 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 100 101 /* Now calculate the permutation and grouping information. */ 102 PetscCall(MatMPIAIJCRL_create_aijcrl(A)); 103 PetscFunctionReturn(0); 104 } 105 106 extern PetscErrorCode MatMult_AIJCRL(Mat,Vec,Vec); 107 extern PetscErrorCode MatDuplicate_AIJCRL(Mat,MatDuplicateOption,Mat*); 108 109 /* MatConvert_MPIAIJ_MPIAIJCRL converts a MPIAIJ matrix into a 110 * MPIAIJCRL matrix. This routine is called by the MatCreate_MPIAIJCRL() 111 * routine, but can also be used to convert an assembled MPIAIJ matrix 112 * into a MPIAIJCRL one. */ 113 114 PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_MPIAIJCRL(Mat A,MatType type,MatReuse reuse,Mat *newmat) 115 { 116 Mat B = *newmat; 117 Mat_AIJCRL *aijcrl; 118 119 PetscFunctionBegin; 120 if (reuse == MAT_INITIAL_MATRIX) { 121 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B)); 122 } 123 124 PetscCall(PetscNewLog(B,&aijcrl)); 125 B->spptr = (void*) aijcrl; 126 127 /* Set function pointers for methods that we inherit from AIJ but override. */ 128 B->ops->duplicate = MatDuplicate_AIJCRL; 129 B->ops->assemblyend = MatAssemblyEnd_MPIAIJCRL; 130 B->ops->destroy = MatDestroy_MPIAIJCRL; 131 B->ops->mult = MatMult_AIJCRL; 132 133 /* If A has already been assembled, compute the permutation. */ 134 if (A->assembled) PetscCall(MatMPIAIJCRL_create_aijcrl(B)); 135 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATMPIAIJCRL)); 136 *newmat = B; 137 PetscFunctionReturn(0); 138 } 139 140 /*@C 141 MatCreateMPIAIJCRL - Creates a sparse matrix of type MPIAIJCRL. 142 This type inherits from AIJ, but stores some additional 143 information that is used to allow better vectorization of 144 the matrix-vector product. At the cost of increased storage, the AIJ formatted 145 matrix can be copied to a format in which pieces of the matrix are 146 stored in ELLPACK format, allowing the vectorized matrix multiply 147 routine to use stride-1 memory accesses. As with the AIJ type, it is 148 important to preallocate matrix storage in order to get good assembly 149 performance. 150 151 Collective 152 153 Input Parameters: 154 + comm - MPI communicator, set to PETSC_COMM_SELF 155 . m - number of rows 156 . n - number of columns 157 . nz - number of nonzeros per row (same for all rows) 158 - nnz - array containing the number of nonzeros in the various rows 159 (possibly different for each row) or NULL 160 161 Output Parameter: 162 . A - the matrix 163 164 Notes: 165 If nnz is given then nz is ignored 166 167 Level: intermediate 168 169 .seealso: `MatCreate()`, `MatCreateMPIAIJPERM()`, `MatSetValues()` 170 @*/ 171 PetscErrorCode MatCreateMPIAIJCRL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],PetscInt onz,const PetscInt onnz[],Mat *A) 172 { 173 PetscFunctionBegin; 174 PetscCall(MatCreate(comm,A)); 175 PetscCall(MatSetSizes(*A,m,n,m,n)); 176 PetscCall(MatSetType(*A,MATMPIAIJCRL)); 177 PetscCall(MatMPIAIJSetPreallocation_MPIAIJ(*A,nz,(PetscInt*)nnz,onz,(PetscInt*)onnz)); 178 PetscFunctionReturn(0); 179 } 180 181 PETSC_EXTERN PetscErrorCode MatCreate_MPIAIJCRL(Mat A) 182 { 183 PetscFunctionBegin; 184 PetscCall(MatSetType(A,MATMPIAIJ)); 185 PetscCall(MatConvert_MPIAIJ_MPIAIJCRL(A,MATMPIAIJCRL,MAT_INPLACE_MATRIX,&A)); 186 PetscFunctionReturn(0); 187 } 188