1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 Mat A; 6 } Mat_Transpose; 7 8 PetscErrorCode MatMult_Transpose(Mat N,Vec x,Vec y) 9 { 10 Mat_Transpose *Na = (Mat_Transpose*)N->data; 11 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 ierr = MatMultTranspose(Na->A,x,y);CHKERRQ(ierr); 15 PetscFunctionReturn(0); 16 } 17 18 PetscErrorCode MatMultAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3) 19 { 20 Mat_Transpose *Na = (Mat_Transpose*)N->data; 21 PetscErrorCode ierr; 22 23 PetscFunctionBegin; 24 ierr = MatMultTransposeAdd(Na->A,v1,v2,v3);CHKERRQ(ierr); 25 PetscFunctionReturn(0); 26 } 27 28 PetscErrorCode MatMultTranspose_Transpose(Mat N,Vec x,Vec y) 29 { 30 Mat_Transpose *Na = (Mat_Transpose*)N->data; 31 PetscErrorCode ierr; 32 33 PetscFunctionBegin; 34 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 38 PetscErrorCode MatMultTransposeAdd_Transpose(Mat N,Vec v1,Vec v2,Vec v3) 39 { 40 Mat_Transpose *Na = (Mat_Transpose*)N->data; 41 PetscErrorCode ierr; 42 43 PetscFunctionBegin; 44 ierr = MatMultAdd(Na->A,v1,v2,v3);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 PetscErrorCode MatDestroy_Transpose(Mat N) 49 { 50 Mat_Transpose *Na = (Mat_Transpose*)N->data; 51 PetscErrorCode ierr; 52 53 PetscFunctionBegin; 54 ierr = MatDestroy(&Na->A);CHKERRQ(ierr); 55 ierr = PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);CHKERRQ(ierr); 56 ierr = PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);CHKERRQ(ierr); 57 ierr = PetscFree(N->data);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m) 62 { 63 Mat_Transpose *Na = (Mat_Transpose*)N->data; 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 if (op == MAT_COPY_VALUES) { 68 ierr = MatTranspose(Na->A,MAT_INITIAL_MATRIX,m);CHKERRQ(ierr); 69 } else if (op == MAT_DO_NOT_COPY_VALUES) { 70 ierr = MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);CHKERRQ(ierr); 71 ierr = MatTranspose(*m,MAT_INPLACE_MATRIX,m);CHKERRQ(ierr); 72 } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type"); 73 PetscFunctionReturn(0); 74 } 75 76 PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l) 77 { 78 Mat_Transpose *Aa = (Mat_Transpose*)A->data; 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 ierr = MatCreateVecs(Aa->A,l,r);CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 PetscErrorCode MatAXPY_Transpose(Mat Y,PetscScalar a,Mat X,MatStructure str) 87 { 88 Mat_Transpose *Ya = (Mat_Transpose*)Y->data; 89 Mat_Transpose *Xa = (Mat_Transpose*)X->data; 90 Mat M = Ya->A; 91 Mat N = Xa->A; 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 ierr = MatAXPY(M,a,N,str);CHKERRQ(ierr); 96 PetscFunctionReturn(0); 97 } 98 99 PetscErrorCode MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool *has) 100 { 101 Mat_Transpose *X = (Mat_Transpose*)mat->data; 102 PetscErrorCode ierr; 103 PetscFunctionBegin; 104 105 *has = PETSC_FALSE; 106 if (op == MATOP_MULT) { 107 ierr = MatHasOperation(X->A,MATOP_MULT_TRANSPOSE,has);CHKERRQ(ierr); 108 } else if (op == MATOP_MULT_TRANSPOSE) { 109 ierr = MatHasOperation(X->A,MATOP_MULT,has);CHKERRQ(ierr); 110 } else if (op == MATOP_MULT_ADD) { 111 ierr = MatHasOperation(X->A,MATOP_MULT_TRANSPOSE_ADD,has);CHKERRQ(ierr); 112 } else if (op == MATOP_MULT_TRANSPOSE_ADD) { 113 ierr = MatHasOperation(X->A,MATOP_MULT_ADD,has);CHKERRQ(ierr); 114 } else if (((void**)mat->ops)[op]) *has = PETSC_TRUE; 115 PetscFunctionReturn(0); 116 } 117 118 /* used by hermitian transpose */ 119 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D) 120 { 121 Mat A,B,C,Ain,Bin,Cin; 122 PetscBool Aistrans,Bistrans,Cistrans; 123 PetscInt Atrans,Btrans,Ctrans; 124 MatProductType ptype; 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 MatCheckProduct(D,1); 129 A = D->product->A; 130 B = D->product->B; 131 C = D->product->C; 132 ierr = PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&Aistrans);CHKERRQ(ierr); 133 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&Bistrans);CHKERRQ(ierr); 134 ierr = PetscObjectTypeCompare((PetscObject)C,MATTRANSPOSEMAT,&Cistrans);CHKERRQ(ierr); 135 if (!Aistrans && !Bistrans && !Cistrans) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"This should not happen"); 136 Atrans = 0; 137 Ain = A; 138 while (Aistrans) { 139 Atrans++; 140 ierr = MatTransposeGetMat(Ain,&Ain);CHKERRQ(ierr); 141 ierr = PetscObjectTypeCompare((PetscObject)Ain,MATTRANSPOSEMAT,&Aistrans);CHKERRQ(ierr); 142 } 143 Btrans = 0; 144 Bin = B; 145 while (Bistrans) { 146 Btrans++; 147 ierr = MatTransposeGetMat(Bin,&Bin);CHKERRQ(ierr); 148 ierr = PetscObjectTypeCompare((PetscObject)Bin,MATTRANSPOSEMAT,&Bistrans);CHKERRQ(ierr); 149 } 150 Ctrans = 0; 151 Cin = C; 152 while (Cistrans) { 153 Ctrans++; 154 ierr = MatTransposeGetMat(Cin,&Cin);CHKERRQ(ierr); 155 ierr = PetscObjectTypeCompare((PetscObject)Cin,MATTRANSPOSEMAT,&Cistrans);CHKERRQ(ierr); 156 } 157 Atrans = Atrans%2; 158 Btrans = Btrans%2; 159 Ctrans = Ctrans%2; 160 ptype = D->product->type; /* same product type by default */ 161 if (Ain->symmetric) Atrans = 0; 162 if (Bin->symmetric) Btrans = 0; 163 if (Cin && Cin->symmetric) Ctrans = 0; 164 165 if (Atrans || Btrans || Ctrans) { 166 ptype = MATPRODUCT_UNSPECIFIED; 167 switch (D->product->type) { 168 case MATPRODUCT_AB: 169 if (Atrans && Btrans) { /* At * Bt we do not have support for this */ 170 /* TODO custom implementation ? */ 171 } else if (Atrans) { /* At * B */ 172 ptype = MATPRODUCT_AtB; 173 } else { /* A * Bt */ 174 ptype = MATPRODUCT_ABt; 175 } 176 break; 177 case MATPRODUCT_AtB: 178 if (Atrans && Btrans) { /* A * Bt */ 179 ptype = MATPRODUCT_ABt; 180 } else if (Atrans) { /* A * B */ 181 ptype = MATPRODUCT_AB; 182 } else { /* At * Bt we do not have support for this */ 183 /* TODO custom implementation ? */ 184 } 185 break; 186 case MATPRODUCT_ABt: 187 if (Atrans && Btrans) { /* At * B */ 188 ptype = MATPRODUCT_AtB; 189 } else if (Atrans) { /* At * Bt we do not have support for this */ 190 /* TODO custom implementation ? */ 191 } else { /* A * B */ 192 ptype = MATPRODUCT_AB; 193 } 194 break; 195 case MATPRODUCT_PtAP: 196 if (Atrans) { /* PtAtP */ 197 /* TODO custom implementation ? */ 198 } else { /* RARt */ 199 ptype = MATPRODUCT_RARt; 200 } 201 break; 202 case MATPRODUCT_RARt: 203 if (Atrans) { /* RAtRt */ 204 /* TODO custom implementation ? */ 205 } else { /* PtAP */ 206 ptype = MATPRODUCT_PtAP; 207 } 208 break; 209 case MATPRODUCT_ABC: 210 /* TODO custom implementation ? */ 211 break; 212 default: SETERRQ1(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"ProductType %s is not supported",MatProductTypes[D->product->type]); 213 } 214 } 215 ierr = MatProductReplaceMats(Ain,Bin,Cin,D);CHKERRQ(ierr); 216 ierr = MatProductSetType(D,ptype);CHKERRQ(ierr); 217 ierr = MatProductSetFromOptions(D);CHKERRQ(ierr); 218 PetscFunctionReturn(0); 219 } 220 221 PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M) 222 { 223 Mat_Transpose *Aa = (Mat_Transpose*)A->data; 224 225 PetscFunctionBegin; 226 *M = Aa->A; 227 PetscFunctionReturn(0); 228 } 229 230 /*@ 231 MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT 232 233 Logically collective on Mat 234 235 Input Parameter: 236 . A - the MATTRANSPOSE matrix 237 238 Output Parameter: 239 . M - the matrix object stored inside A 240 241 Level: intermediate 242 243 .seealso: MatCreateTranspose() 244 245 @*/ 246 PetscErrorCode MatTransposeGetMat(Mat A,Mat *M) 247 { 248 PetscErrorCode ierr; 249 250 PetscFunctionBegin; 251 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 252 PetscValidType(A,1); 253 PetscValidPointer(M,2); 254 ierr = PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 258 /*@ 259 MatCreateTranspose - Creates a new matrix object that behaves like A' 260 261 Collective on Mat 262 263 Input Parameter: 264 . A - the (possibly rectangular) matrix 265 266 Output Parameter: 267 . N - the matrix that represents A' 268 269 Level: intermediate 270 271 Notes: 272 The transpose A' is NOT actually formed! Rather the new matrix 273 object performs the matrix-vector product by using the MatMultTranspose() on 274 the original matrix 275 276 .seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate() 277 278 @*/ 279 PetscErrorCode MatCreateTranspose(Mat A,Mat *N) 280 { 281 PetscErrorCode ierr; 282 PetscInt m,n; 283 Mat_Transpose *Na; 284 VecType vtype; 285 286 PetscFunctionBegin; 287 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 288 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr); 289 ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 290 ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr); 291 ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr); 292 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr); 293 294 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 295 (*N)->data = (void*) Na; 296 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 297 Na->A = A; 298 299 (*N)->ops->destroy = MatDestroy_Transpose; 300 (*N)->ops->mult = MatMult_Transpose; 301 (*N)->ops->multadd = MatMultAdd_Transpose; 302 (*N)->ops->multtranspose = MatMultTranspose_Transpose; 303 (*N)->ops->multtransposeadd = MatMultTransposeAdd_Transpose; 304 (*N)->ops->duplicate = MatDuplicate_Transpose; 305 (*N)->ops->getvecs = MatCreateVecs_Transpose; 306 (*N)->ops->axpy = MatAXPY_Transpose; 307 (*N)->ops->hasoperation = MatHasOperation_Transpose; 308 (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose; 309 (*N)->assembled = PETSC_TRUE; 310 311 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);CHKERRQ(ierr); 312 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);CHKERRQ(ierr); 313 ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 314 ierr = MatGetVecType(A,&vtype);CHKERRQ(ierr); 315 ierr = MatSetVecType(*N,vtype);CHKERRQ(ierr); 316 ierr = MatSetUp(*N);CHKERRQ(ierr); 317 PetscFunctionReturn(0); 318 } 319