1 /* 2 Defines basic operations for the MATSEQAIJSELL matrix class. 3 This class is derived from the MATAIJCLASS, but maintains a "shadow" copy 4 of the matrix stored in MATSEQSELL format, which is used as appropriate for 5 performing operations for which this format is more suitable. 6 */ 7 8 #include <../src/mat/impls/aij/seq/aij.h> 9 10 typedef struct { 11 Mat S; /* The SELL formatted "shadow" matrix. */ 12 PetscBool eager_shadow; 13 PetscObjectState state; /* State of the matrix when shadow matrix was last constructed. */ 14 } Mat_SeqAIJSELL; 15 16 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJSELL_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat) 17 { 18 /* This routine is only called to convert a MATAIJSELL to its base PETSc type, */ 19 /* so we will ignore 'MatType type'. */ 20 PetscErrorCode ierr; 21 Mat B = *newmat; 22 Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL*) A->spptr; 23 24 PetscFunctionBegin; 25 if (reuse == MAT_INITIAL_MATRIX) { 26 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 27 } 28 29 /* Reset the original function pointers. */ 30 B->ops->duplicate = MatDuplicate_SeqAIJ; 31 B->ops->assemblyend = MatAssemblyEnd_SeqAIJ; 32 B->ops->destroy = MatDestroy_SeqAIJ; 33 B->ops->mult = MatMult_SeqAIJ; 34 B->ops->multtranspose = MatMultTranspose_SeqAIJ; 35 B->ops->multadd = MatMultAdd_SeqAIJ; 36 B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJ; 37 38 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaijsell_seqaij_C",NULL);CHKERRQ(ierr); 39 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaijsell_C",NULL);CHKERRQ(ierr); 40 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaijsell_C",NULL);CHKERRQ(ierr); 41 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaijsell_C",NULL);CHKERRQ(ierr); 42 43 if (reuse == MAT_INITIAL_MATRIX) aijsell = (Mat_SeqAIJSELL*)B->spptr; 44 45 /* Clean up the Mat_SeqAIJSELL data structure. */ 46 if(aijsell->S) { 47 ierr = MatDestroy(&aijsell->S);CHKERRQ(ierr); 48 } 49 ierr = PetscFree(B->spptr);CHKERRQ(ierr); 50 51 /* Change the type of B to MATSEQAIJ. */ 52 ierr = PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ);CHKERRQ(ierr); 53 54 *newmat = B; 55 PetscFunctionReturn(0); 56 } 57 58 PetscErrorCode MatDestroy_SeqAIJSELL(Mat A) 59 { 60 PetscErrorCode ierr; 61 Mat_SeqAIJSELL *aijsell = (Mat_SeqAIJSELL*) A->spptr; 62 63 PetscFunctionBegin; 64 65 /* If MatHeaderMerge() was used, then this SeqAIJSELL matrix will not have an 66 * spptr pointer. */ 67 if (aijsell) { 68 /* Clean up everything in the Mat_SeqAIJSELL data structure, then free A->spptr. */ 69 if (aijsell->S) { 70 ierr = MatDestroy(&aijsell->S);CHKERRQ(ierr); 71 } 72 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 73 } 74 75 /* Change the type of A back to SEQAIJ and use MatDestroy_SeqAIJ() 76 * to destroy everything that remains. */ 77 ierr = PetscObjectChangeTypeName((PetscObject)A, MATSEQAIJ);CHKERRQ(ierr); 78 /* Note that I don't call MatSetType(). I believe this is because that 79 * is only to be called when *building* a matrix. I could be wrong, but 80 * that is how things work for the SuperLU matrix class. */ 81 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 PetscErrorCode MatDuplicate_SeqAIJSELL(Mat A, MatDuplicateOption op, Mat *M) 86 { 87 PetscErrorCode ierr; 88 Mat_SeqAIJSELL *aijsell; 89 Mat_SeqAIJSELL *aijsell_dest; 90 91 PetscFunctionBegin; 92 ierr = MatDuplicate_SeqAIJ(A,op,M);CHKERRQ(ierr); 93 aijsell = (Mat_SeqAIJSELL*) A->spptr; 94 aijsell_dest = (Mat_SeqAIJSELL*) (*M)->spptr; 95 ierr = PetscMemcpy(aijsell_dest,aijsell,sizeof(Mat_SeqAIJSELL));CHKERRQ(ierr); 96 /* We don't duplicate the shadow matrix -- that will be constructed as needed. */ 97 aijsell_dest->S = NULL; 98 /* TODO: Have the shadow matrix be built now if eager_shadow is set to true. */ 99 PetscFunctionReturn(0); 100 } 101 102 PetscErrorCode MatAssemblyEnd_SeqAIJSELL(Mat A, MatAssemblyType mode) 103 { 104 PetscErrorCode ierr; 105 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 106 Mat_SeqAIJSELL *aijsell; 107 108 PetscFunctionBegin; 109 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 110 111 /* I disable the use of the inode routines so that the AIJSELL ones will be 112 * used instead, but I wonder if it might make sense (and is feasible) to 113 * use some of them. */ 114 a->inode.use = PETSC_FALSE; 115 116 /* Since a MATSEQAIJSELL matrix is really just a MATSEQAIJ with some 117 * extra information and some different methods, call the AssemblyEnd 118 * routine for a MATSEQAIJ. 119 * I'm not sure if this is the best way to do this, but it avoids 120 * a lot of code duplication. */ 121 122 ierr = MatAssemblyEnd_SeqAIJ(A, mode);CHKERRQ(ierr); 123 124 /* If the user has requested "eager" shadowing, create the SELL shadow matrix (if needed; the function checks). 125 * (The default is to take a "lazy" approach, deferring this until something like MatMult() is called.) */ 126 aijsell = (Mat_SeqAIJSELL*) A->spptr; 127 if (aijsell->eager_shadow) { 128 /* TODO: Construct the shadow matrix here. */ 129 } 130 131 PetscFunctionReturn(0); 132 } 133 134 /* MatConvert_SeqAIJ_SeqAIJSELL converts a SeqAIJ matrix into a 135 * SeqAIJSELL matrix. This routine is called by the MatCreate_SeqAIJSELL() 136 * routine, but can also be used to convert an assembled SeqAIJ matrix 137 * into a SeqAIJSELL one. */ 138 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJSELL(Mat A,MatType type,MatReuse reuse,Mat *newmat) 139 { 140 PetscErrorCode ierr; 141 Mat B = *newmat; 142 Mat_SeqAIJ *b; 143 Mat_SeqAIJSELL *aijsell; 144 PetscBool set; 145 PetscBool sametype; 146 147 PetscFunctionBegin; 148 if (reuse == MAT_INITIAL_MATRIX) { 149 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 150 } 151 152 ierr = PetscObjectTypeCompare((PetscObject)A,type,&sametype);CHKERRQ(ierr); 153 if (sametype) PetscFunctionReturn(0); 154 155 ierr = PetscNewLog(B,&aijsell);CHKERRQ(ierr); 156 b = (Mat_SeqAIJ*) B->data; 157 B->spptr = (void*) aijsell; 158 159 /* Disable use of the inode routines so that the AIJSELL ones will be used instead. 160 * This happens in MatAssemblyEnd_SeqAIJSELL as well, but the assembly end may not be called, so set it here, too. 161 * As noted elsewhere, I wonder if it might make sense and be feasible to use some of the inode routines. */ 162 b->inode.use = PETSC_FALSE; 163 164 /* Set function pointers for methods that we inherit from AIJ but override. 165 * We also parse some command line options below, since those determine some of the methods we point to. */ 166 B->ops->duplicate = MatDuplicate_SeqAIJSELL; 167 B->ops->assemblyend = MatAssemblyEnd_SeqAIJSELL; 168 B->ops->destroy = MatDestroy_SeqAIJSELL; 169 170 aijsell->S = NULL; 171 aijsell->eager_shadow = PETSC_FALSE; 172 173 /* Parse command line options. */ 174 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"AIJSELL Options","Mat");CHKERRQ(ierr); 175 ierr = PetscOptionsBool("-mat_aijsell_eager_shadow","Eager Shadowing","None",(PetscBool)aijsell->eager_shadow,(PetscBool*)&aijsell->eager_shadow,&set);CHKERRQ(ierr); 176 ierr = PetscOptionsEnd();CHKERRQ(ierr); 177 178 /* If A has already been assembled and eager shadowing is specified, build the shadow matrix. */ 179 if (A->assembled) { 180 /* TODO: Have the shadow matrix for B be built now if eager_shadow is set to true. */ 181 } 182 183 /* Commenting these out for now, as these functions are not yet implemented. 184 B->ops->mult = MatMult_SeqAIJSELL_SpMV2; 185 B->ops->multtranspose = MatMultTranspose_SeqAIJSELL_SpMV2; 186 B->ops->multadd = MatMultAdd_SeqAIJSELL_SpMV2; 187 B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJSELL_SpMV2; 188 */ 189 190 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaijsell_seqaij_C",MatConvert_SeqAIJSELL_SeqAIJ);CHKERRQ(ierr); 191 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaijsell_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 192 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaijsell_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 193 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaijsell_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 194 195 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJSELL);CHKERRQ(ierr); 196 *newmat = B; 197 PetscFunctionReturn(0); 198 } 199 200 /*@C 201 MatCreateSeqAIJSELL - Creates a sparse matrix of type SEQAIJSELL. 202 This type inherits from AIJ and is largely identical, but keeps a "shadow" 203 copy of the matrix in SEQSELL format, which is used when performing 204 operations for which this format is more suitable. 205 206 Collective on MPI_Comm 207 208 Input Parameters: 209 + comm - MPI communicator, set to PETSC_COMM_SELF 210 . m - number of rows 211 . n - number of columns 212 . nz - number of nonzeros per row (same for all rows) 213 - nnz - array containing the number of nonzeros in the various rows 214 (possibly different for each row) or NULL 215 216 Output Parameter: 217 . A - the matrix 218 219 Options Database Keys: 220 . -mat_aijsell_eager_shadow - Construct shadow matrix upon matrix assembly; default is to take a "lazy" approach, performing this step the first time the matrix is applied 221 222 Notes: 223 If nnz is given then nz is ignored 224 225 Level: intermediate 226 227 .keywords: matrix, sparse, parallel 228 229 .seealso: MatCreate(), MatCreateMPIAIJSELL(), MatSetValues() 230 @*/ 231 PetscErrorCode MatCreateSeqAIJSELL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 232 { 233 PetscErrorCode ierr; 234 235 PetscFunctionBegin; 236 ierr = MatCreate(comm,A);CHKERRQ(ierr); 237 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 238 ierr = MatSetType(*A,MATSEQAIJSELL);CHKERRQ(ierr); 239 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 243 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJSELL(Mat A) 244 { 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 249 ierr = MatConvert_SeqAIJ_SeqAIJSELL(A,MATSEQAIJSELL,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252