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 PetscErrorCode ierr; 95 PetscTruth flg; 96 97 PetscFunctionBegin; 98 ierr = PetscOptionsHasName(Y->prefix,"-mat_composite_merge",&flg);CHKERRQ(ierr); 99 if (flg) { 100 ierr = MatCompositeMerge(Y);CHKERRQ(ierr); 101 } 102 PetscFunctionReturn(0); 103 } 104 105 106 static struct _MatOps MatOps_Values = {0, 107 0, 108 0, 109 MatMult_Composite, 110 0, 111 /* 5*/ MatMultTranspose_Composite, 112 0, 113 0, 114 0, 115 0, 116 /*10*/ 0, 117 0, 118 0, 119 0, 120 0, 121 /*15*/ 0, 122 0, 123 MatGetDiagonal_Composite, 124 0, 125 0, 126 /*20*/ 0, 127 MatAssemblyEnd_Composite, 128 0, 129 0, 130 0, 131 /*25*/ 0, 132 0, 133 0, 134 0, 135 0, 136 /*30*/ 0, 137 0, 138 0, 139 0, 140 0, 141 /*35*/ 0, 142 0, 143 0, 144 0, 145 0, 146 /*40*/ 0, 147 0, 148 0, 149 0, 150 0, 151 /*45*/ 0, 152 0, 153 0, 154 0, 155 0, 156 /*50*/ 0, 157 0, 158 0, 159 0, 160 0, 161 /*55*/ 0, 162 0, 163 0, 164 0, 165 0, 166 /*60*/ 0, 167 MatDestroy_Composite, 168 0, 169 0, 170 0, 171 /*65*/ 0, 172 0, 173 0, 174 0, 175 0, 176 /*70*/ 0, 177 0, 178 0, 179 0, 180 0, 181 /*75*/ 0, 182 0, 183 0, 184 0, 185 0, 186 /*80*/ 0, 187 0, 188 0, 189 0, 190 0, 191 /*85*/ 0, 192 0, 193 0, 194 0, 195 0, 196 /*90*/ 0, 197 0, 198 0, 199 0, 200 0, 201 /*95*/ 0, 202 0, 203 0, 204 0}; 205 206 /*MC 207 MATCOMPOSITE - A matrix defined by the sum of one or more matrices (all matrices are of same size and parallel layout) 208 209 Level: advanced 210 211 .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge() 212 M*/ 213 214 EXTERN_C_BEGIN 215 #undef __FUNCT__ 216 #define __FUNCT__ "MatCreate_Composite" 217 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Composite(Mat A) 218 { 219 Mat_Composite *b; 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 224 225 ierr = PetscNew(Mat_Composite,&b);CHKERRQ(ierr); 226 A->data = (void*)b; 227 228 ierr = PetscMapSetBlockSize(&A->rmap,1);CHKERRQ(ierr); 229 ierr = PetscMapSetBlockSize(&A->cmap,1);CHKERRQ(ierr); 230 ierr = PetscMapSetUp(&A->rmap);CHKERRQ(ierr); 231 ierr = PetscMapSetUp(&A->cmap);CHKERRQ(ierr); 232 233 A->assembled = PETSC_TRUE; 234 A->preallocated = PETSC_FALSE; 235 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr); 236 PetscFunctionReturn(0); 237 } 238 EXTERN_C_END 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "MatCreateComposite" 242 /*@C 243 MatCreateComposite - Creates a matrix as the sum of zero or more matrices 244 245 Collective on MPI_Comm 246 247 Input Parameters: 248 + comm - MPI communicator 249 . nmat - number of matrices to put in 250 - mats - the matrices 251 252 Output Parameter: 253 . A - the matrix 254 255 Level: advanced 256 257 Notes: 258 Alternative construction 259 $ MatCreate(comm,&mat); 260 $ MatSetSizes(mat,m,n,M,N); 261 $ MatSetType(mat,MATCOMPOSITE); 262 $ MatCompositeAddMat(mat,mats[0]); 263 $ .... 264 $ MatCompositeAddMat(mat,mats[nmat-1]); 265 $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 266 $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); 267 268 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge() 269 270 @*/ 271 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat) 272 { 273 PetscErrorCode ierr; 274 PetscInt m,n,M,N,i; 275 276 PetscFunctionBegin; 277 if (nmat < 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix"); 278 PetscValidPointer(mat,3); 279 280 ierr = MatGetLocalSize(mats[0],&m,&n);CHKERRQ(ierr); 281 ierr = MatGetSize(mats[0],&M,&N);CHKERRQ(ierr); 282 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 283 ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr); 284 ierr = MatSetType(*mat,MATCOMPOSITE);CHKERRQ(ierr); 285 for (i=0; i<nmat; i++) { 286 ierr = MatCompositeAddMat(*mat,mats[i]);CHKERRQ(ierr); 287 } 288 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 289 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 #undef __FUNCT__ 294 #define __FUNCT__ "MatCompositeAddMat" 295 /*@ 296 MatCompositeAddMat - add another matrix to a composite matrix 297 298 Collective on Mat 299 300 Input Parameters: 301 + mat - the composite matrix 302 - smat - the partial matrix 303 304 Level: advanced 305 306 .seealso: MatCreateComposite() 307 @*/ 308 PetscErrorCode PETSCMAT_DLLEXPORT MatCompositeAddMat(Mat mat,Mat smat) 309 { 310 Mat_Composite *shell = (Mat_Composite*)mat->data; 311 PetscErrorCode ierr; 312 Mat_CompositeLink ilink,next; 313 314 PetscFunctionBegin; 315 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 316 PetscValidHeaderSpecific(smat,MAT_COOKIE,2); 317 ierr = PetscNew(struct _Mat_CompositeLink,&ilink);CHKERRQ(ierr); 318 ilink->next = 0; 319 ierr = PetscObjectReference((PetscObject)smat);CHKERRQ(ierr); 320 ilink->mat = smat; 321 322 next = shell->head; 323 if (!next) { 324 shell->head = ilink; 325 } else { 326 while (next->next) { 327 next = next->next; 328 } 329 next->next = ilink; 330 } 331 PetscFunctionReturn(0); 332 } 333 334 #undef __FUNCT__ 335 #define __FUNCT__ "MatCompositeMerge" 336 /*@C 337 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 338 by summing all the matrices inside the composite matrix. 339 340 Collective on MPI_Comm 341 342 Input Parameters: 343 . mat - the composite matrix 344 345 346 Options Database: 347 . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked) 348 349 Level: advanced 350 351 Notes: 352 The MatType of the resulting matrix will be the same as the MatType of the FIRST 353 matrix in the composite matrix. 354 355 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE 356 357 @*/ 358 PetscErrorCode PETSCMAT_DLLEXPORT MatCompositeMerge(Mat mat) 359 { 360 Mat_Composite *shell = (Mat_Composite*)mat->data; 361 Mat_CompositeLink next = shell->head; 362 PetscErrorCode ierr; 363 Mat tmat; 364 365 PetscFunctionBegin; 366 if (!next) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 367 368 PetscFunctionBegin; 369 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 370 while ((next = next->next)) { 371 ierr = MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 372 } 373 ierr = MatDestroy_Composite(mat);CHKERRQ(ierr); 374 375 PetscFunctionReturn(0); 376 } 377