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) { 135 PetscCall(MatMPIAIJCRL_create_aijcrl(B)); 136 } 137 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATMPIAIJCRL)); 138 *newmat = B; 139 PetscFunctionReturn(0); 140 } 141 142 /*@C 143 MatCreateMPIAIJCRL - Creates a sparse matrix of type MPIAIJCRL. 144 This type inherits from AIJ, but stores some additional 145 information that is used to allow better vectorization of 146 the matrix-vector product. At the cost of increased storage, the AIJ formatted 147 matrix can be copied to a format in which pieces of the matrix are 148 stored in ELLPACK format, allowing the vectorized matrix multiply 149 routine to use stride-1 memory accesses. As with the AIJ type, it is 150 important to preallocate matrix storage in order to get good assembly 151 performance. 152 153 Collective 154 155 Input Parameters: 156 + comm - MPI communicator, set to PETSC_COMM_SELF 157 . m - number of rows 158 . n - number of columns 159 . nz - number of nonzeros per row (same for all rows) 160 - nnz - array containing the number of nonzeros in the various rows 161 (possibly different for each row) or NULL 162 163 Output Parameter: 164 . A - the matrix 165 166 Notes: 167 If nnz is given then nz is ignored 168 169 Level: intermediate 170 171 .seealso: MatCreate(), MatCreateMPIAIJPERM(), MatSetValues() 172 @*/ 173 PetscErrorCode MatCreateMPIAIJCRL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],PetscInt onz,const PetscInt onnz[],Mat *A) 174 { 175 PetscFunctionBegin; 176 PetscCall(MatCreate(comm,A)); 177 PetscCall(MatSetSizes(*A,m,n,m,n)); 178 PetscCall(MatSetType(*A,MATMPIAIJCRL)); 179 PetscCall(MatMPIAIJSetPreallocation_MPIAIJ(*A,nz,(PetscInt*)nnz,onz,(PetscInt*)onnz)); 180 PetscFunctionReturn(0); 181 } 182 183 PETSC_EXTERN PetscErrorCode MatCreate_MPIAIJCRL(Mat A) 184 { 185 PetscFunctionBegin; 186 PetscCall(MatSetType(A,MATMPIAIJ)); 187 PetscCall(MatConvert_MPIAIJ_MPIAIJCRL(A,MATMPIAIJCRL,MAT_INPLACE_MATRIX,&A)); 188 PetscFunctionReturn(0); 189 } 190