1 #define PETSCMAT_DLL 2 3 #include "include/private/matimpl.h" /*I "petscmat.h" I*/ 4 5 typedef struct _Mat_CompositeLink *Mat_CompositeLink; 6 struct _Mat_CompositeLink { 7 Mat mat; 8 Mat_CompositeLink next; 9 }; 10 11 typedef struct { 12 Mat_CompositeLink head; 13 Vec work; 14 } Mat_Composite; 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "MatDestroy_Composite" 18 PetscErrorCode MatDestroy_Composite(Mat mat) 19 { 20 PetscErrorCode ierr; 21 Mat_Composite *shell = (Mat_Composite*)mat->data; 22 Mat_CompositeLink next = shell->head; 23 24 PetscFunctionBegin; 25 while (next) { 26 ierr = MatDestroy(next->mat);CHKERRQ(ierr); 27 next = next->next; 28 } 29 if (shell->work) {ierr = VecDestroy(shell->work);CHKERRQ(ierr);} 30 ierr = PetscFree(shell);CHKERRQ(ierr); 31 mat->data = 0; 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "MatMult_Composite" 37 PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y) 38 { 39 Mat_Composite *shell = (Mat_Composite*)A->data; 40 Mat_CompositeLink next = shell->head; 41 PetscErrorCode ierr; 42 43 PetscFunctionBegin; 44 if (!next) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 45 ierr = MatMult(next->mat,x,y);CHKERRQ(ierr); 46 while ((next = next->next)) { 47 ierr = MatMultAdd(next->mat,x,y,y);CHKERRQ(ierr); 48 } 49 PetscFunctionReturn(0); 50 } 51 52 #undef __FUNCT__ 53 #define __FUNCT__ "MatMultTranspose_Composite" 54 PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y) 55 { 56 Mat_Composite *shell = (Mat_Composite*)A->data; 57 Mat_CompositeLink next = shell->head; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 if (!next) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 62 ierr = MatMultTranspose(next->mat,x,y);CHKERRQ(ierr); 63 while ((next = next->next)) { 64 ierr = MatMultTransposeAdd(next->mat,x,y,y);CHKERRQ(ierr); 65 } 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "MatGetDiagonal_Composite" 71 PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v) 72 { 73 Mat_Composite *shell = (Mat_Composite*)A->data; 74 Mat_CompositeLink next = shell->head; 75 PetscErrorCode ierr; 76 77 PetscFunctionBegin; 78 if (!next) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 79 ierr = MatGetDiagonal(next->mat,v);CHKERRQ(ierr); 80 if (next->next && !shell->work) { 81 ierr = VecDuplicate(v,&shell->work);CHKERRQ(ierr); 82 } 83 while ((next = next->next)) { 84 ierr = MatGetDiagonal(next->mat,shell->work);CHKERRQ(ierr); 85 ierr = VecAXPY(v,1.0,shell->work);CHKERRQ(ierr); 86 } 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "MatAssemblyEnd_Composite" 92 PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t) 93 { 94 PetscFunctionBegin; 95 PetscFunctionReturn(0); 96 } 97 98 99 static struct _MatOps MatOps_Values = {0, 100 0, 101 0, 102 MatMult_Composite, 103 0, 104 /* 5*/ MatMultTranspose_Composite, 105 0, 106 0, 107 0, 108 0, 109 /*10*/ 0, 110 0, 111 0, 112 0, 113 0, 114 /*15*/ 0, 115 0, 116 MatGetDiagonal_Composite, 117 0, 118 0, 119 /*20*/ 0, 120 MatAssemblyEnd_Composite, 121 0, 122 0, 123 0, 124 /*25*/ 0, 125 0, 126 0, 127 0, 128 0, 129 /*30*/ 0, 130 0, 131 0, 132 0, 133 0, 134 /*35*/ 0, 135 0, 136 0, 137 0, 138 0, 139 /*40*/ 0, 140 0, 141 0, 142 0, 143 0, 144 /*45*/ 0, 145 0, 146 0, 147 0, 148 0, 149 /*50*/ 0, 150 0, 151 0, 152 0, 153 0, 154 /*55*/ 0, 155 0, 156 0, 157 0, 158 0, 159 /*60*/ 0, 160 MatDestroy_Composite, 161 0, 162 0, 163 0, 164 /*65*/ 0, 165 0, 166 0, 167 0, 168 0, 169 /*70*/ 0, 170 0, 171 0, 172 0, 173 0, 174 /*75*/ 0, 175 0, 176 0, 177 0, 178 0, 179 /*80*/ 0, 180 0, 181 0, 182 0, 183 0, 184 /*85*/ 0, 185 0, 186 0, 187 0, 188 0, 189 /*90*/ 0, 190 0, 191 0, 192 0, 193 0, 194 /*95*/ 0, 195 0, 196 0, 197 0}; 198 199 /*MC 200 MATCOMPOSITE - A matrix defined by the sum of one or more matrices (all matrices are of same size and parallel layout) 201 202 Level: advanced 203 204 .seealso: MatCreateComposite 205 M*/ 206 207 EXTERN_C_BEGIN 208 #undef __FUNCT__ 209 #define __FUNCT__ "MatCreate_Composite" 210 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Composite(Mat A) 211 { 212 Mat_Composite *b; 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 217 218 ierr = PetscNew(Mat_Composite,&b);CHKERRQ(ierr); 219 A->data = (void*)b; 220 221 ierr = PetscMapInitialize(A->comm,&A->rmap);CHKERRQ(ierr); 222 ierr = PetscMapInitialize(A->comm,&A->cmap);CHKERRQ(ierr); 223 224 A->assembled = PETSC_TRUE; 225 A->preallocated = PETSC_FALSE; 226 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 EXTERN_C_END 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "MatCreateComposite" 233 /*@C 234 MatCreateComposite - Creates a matrix as the sum of zero or more matrices 235 236 Collective on MPI_Comm 237 238 Input Parameters: 239 + comm - MPI communicator 240 . nmat - number of matrices to put in 241 - mats - the matrices 242 243 Output Parameter: 244 . A - the matrix 245 246 Level: advanced 247 248 Notes: 249 Alternative construction 250 $ MatCreate(comm,&mat); 251 $ MatSetSizes(mat,m,n,M,N); 252 $ MatSetType(mat,MATCOMPOSITE); 253 $ MatCompositeAddMat(mat,mats[0]); 254 $ .... 255 $ MatCompositeAddMat(mat,mats[nmat-1]); 256 257 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat() 258 259 @*/ 260 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat) 261 { 262 PetscErrorCode ierr; 263 PetscInt m,n,M,N,i; 264 265 PetscFunctionBegin; 266 if (nmat < 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix"); 267 PetscValidPointer(mat,3); 268 269 ierr = MatGetLocalSize(mats[0],&m,&n);CHKERRQ(ierr); 270 ierr = MatGetSize(mats[0],&M,&N);CHKERRQ(ierr); 271 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 272 ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr); 273 ierr = MatSetType(*mat,MATCOMPOSITE);CHKERRQ(ierr); 274 for (i=0; i<nmat; i++) { 275 ierr = MatCompositeAddMat(*mat,mats[i]);CHKERRQ(ierr); 276 } 277 PetscFunctionReturn(0); 278 } 279 280 #undef __FUNCT__ 281 #define __FUNCT__ "MatCompositeAddMat" 282 /*@ 283 MatCompositeAddMat - add another matrix to a composite matrix 284 285 Collective on Mat 286 287 Input Parameters: 288 + mat - the composite matrix 289 - smat - the partial matrix 290 291 Level: advanced 292 293 .seealso: MatCreateComposite() 294 @*/ 295 PetscErrorCode PETSCMAT_DLLEXPORT MatCompositeAddMat(Mat mat,Mat smat) 296 { 297 Mat_Composite *shell = (Mat_Composite*)mat->data; 298 PetscErrorCode ierr; 299 Mat_CompositeLink ilink,next; 300 301 PetscFunctionBegin; 302 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 303 PetscValidHeaderSpecific(smat,MAT_COOKIE,2); 304 ierr = PetscNew(struct _Mat_CompositeLink,&ilink);CHKERRQ(ierr); 305 ilink->next = 0; 306 ilink->mat = smat; 307 ierr = PetscObjectReference((PetscObject)smat);CHKERRQ(ierr); 308 309 next = shell->head; 310 if (!next) { 311 shell->head = ilink; 312 } else { 313 while (next->next) { 314 next = next->next; 315 } 316 next->next = ilink; 317 } 318 PetscFunctionReturn(0); 319 } 320 321