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 = PetscMapInitialize(A->comm,&A->rmap);CHKERRQ(ierr); 229 ierr = PetscMapInitialize(A->comm,&A->cmap);CHKERRQ(ierr); 230 231 A->assembled = PETSC_TRUE; 232 A->preallocated = PETSC_FALSE; 233 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr); 234 PetscFunctionReturn(0); 235 } 236 EXTERN_C_END 237 238 #undef __FUNCT__ 239 #define __FUNCT__ "MatCreateComposite" 240 /*@C 241 MatCreateComposite - Creates a matrix as the sum of zero or more matrices 242 243 Collective on MPI_Comm 244 245 Input Parameters: 246 + comm - MPI communicator 247 . nmat - number of matrices to put in 248 - mats - the matrices 249 250 Output Parameter: 251 . A - the matrix 252 253 Level: advanced 254 255 Notes: 256 Alternative construction 257 $ MatCreate(comm,&mat); 258 $ MatSetSizes(mat,m,n,M,N); 259 $ MatSetType(mat,MATCOMPOSITE); 260 $ MatCompositeAddMat(mat,mats[0]); 261 $ .... 262 $ MatCompositeAddMat(mat,mats[nmat-1]); 263 $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 264 $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); 265 266 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge() 267 268 @*/ 269 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat) 270 { 271 PetscErrorCode ierr; 272 PetscInt m,n,M,N,i; 273 274 PetscFunctionBegin; 275 if (nmat < 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix"); 276 PetscValidPointer(mat,3); 277 278 ierr = MatGetLocalSize(mats[0],&m,&n);CHKERRQ(ierr); 279 ierr = MatGetSize(mats[0],&M,&N);CHKERRQ(ierr); 280 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 281 ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr); 282 ierr = MatSetType(*mat,MATCOMPOSITE);CHKERRQ(ierr); 283 for (i=0; i<nmat; i++) { 284 ierr = MatCompositeAddMat(*mat,mats[i]);CHKERRQ(ierr); 285 } 286 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 287 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 } 290 291 #undef __FUNCT__ 292 #define __FUNCT__ "MatCompositeAddMat" 293 /*@ 294 MatCompositeAddMat - add another matrix to a composite matrix 295 296 Collective on Mat 297 298 Input Parameters: 299 + mat - the composite matrix 300 - smat - the partial matrix 301 302 Level: advanced 303 304 .seealso: MatCreateComposite() 305 @*/ 306 PetscErrorCode PETSCMAT_DLLEXPORT MatCompositeAddMat(Mat mat,Mat smat) 307 { 308 Mat_Composite *shell = (Mat_Composite*)mat->data; 309 PetscErrorCode ierr; 310 Mat_CompositeLink ilink,next; 311 312 PetscFunctionBegin; 313 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 314 PetscValidHeaderSpecific(smat,MAT_COOKIE,2); 315 ierr = PetscNew(struct _Mat_CompositeLink,&ilink);CHKERRQ(ierr); 316 ilink->next = 0; 317 ierr = PetscObjectReference((PetscObject)smat);CHKERRQ(ierr); 318 ilink->mat = smat; 319 320 next = shell->head; 321 if (!next) { 322 shell->head = ilink; 323 } else { 324 while (next->next) { 325 next = next->next; 326 } 327 next->next = ilink; 328 } 329 PetscFunctionReturn(0); 330 } 331 332 #undef __FUNCT__ 333 #define __FUNCT__ "MatCompositeMerge" 334 /*@C 335 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 336 by summing all the matrices inside the composite matrix. 337 338 Collective on MPI_Comm 339 340 Input Parameters: 341 . mat - the composite matrix 342 343 344 Options Database: 345 . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked) 346 347 Level: advanced 348 349 Notes: 350 The MatType of the resulting matrix will be the same as the MatType of the FIRST 351 matrix in the composite matrix. 352 353 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE 354 355 @*/ 356 PetscErrorCode PETSCMAT_DLLEXPORT MatCompositeMerge(Mat mat) 357 { 358 Mat_Composite *shell = (Mat_Composite*)mat->data; 359 Mat_CompositeLink next = shell->head; 360 PetscErrorCode ierr; 361 Mat tmat; 362 363 PetscFunctionBegin; 364 if (!next) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 365 366 PetscFunctionBegin; 367 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 368 while ((next = next->next)) { 369 ierr = MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 370 } 371 ierr = MatDestroy_Composite(mat);CHKERRQ(ierr); 372 373 PetscFunctionReturn(0); 374 } 375