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 = PetscFree(N->data);CHKERRQ(ierr); 57 PetscFunctionReturn(0); 58 } 59 60 PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m) 61 { 62 Mat_Transpose *Na = (Mat_Transpose*)N->data; 63 PetscErrorCode ierr; 64 65 PetscFunctionBegin; 66 if (op == MAT_COPY_VALUES) { 67 ierr = MatTranspose(Na->A,MAT_INITIAL_MATRIX,m);CHKERRQ(ierr); 68 } else if (op == MAT_DO_NOT_COPY_VALUES) { 69 ierr = MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);CHKERRQ(ierr); 70 ierr = MatTranspose(*m,MAT_INPLACE_MATRIX,m);CHKERRQ(ierr); 71 } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type"); 72 PetscFunctionReturn(0); 73 } 74 75 PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l) 76 { 77 Mat_Transpose *Aa = (Mat_Transpose*)A->data; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 ierr = MatCreateVecs(Aa->A,l,r);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 PetscErrorCode MatAXPY_Transpose(Mat Y,PetscScalar a,Mat X,MatStructure str) 86 { 87 Mat_Transpose *Ya = (Mat_Transpose*)Y->data; 88 Mat_Transpose *Xa = (Mat_Transpose*)X->data; 89 Mat M = Ya->A; 90 Mat N = Xa->A; 91 PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 ierr = MatAXPY(M,a,N,str);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 98 PetscErrorCode MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool *has) 99 { 100 Mat_Transpose *X = (Mat_Transpose*)mat->data; 101 PetscErrorCode ierr; 102 PetscFunctionBegin; 103 104 *has = PETSC_FALSE; 105 if (op == MATOP_MULT) { 106 ierr = MatHasOperation(X->A,MATOP_MULT_TRANSPOSE,has);CHKERRQ(ierr); 107 } else if (op == MATOP_MULT_TRANSPOSE) { 108 ierr = MatHasOperation(X->A,MATOP_MULT,has);CHKERRQ(ierr); 109 } else if (op == MATOP_MULT_ADD) { 110 ierr = MatHasOperation(X->A,MATOP_MULT_TRANSPOSE_ADD,has);CHKERRQ(ierr); 111 } else if (op == MATOP_MULT_TRANSPOSE_ADD) { 112 ierr = MatHasOperation(X->A,MATOP_MULT_ADD,has);CHKERRQ(ierr); 113 } else if (((void**)mat->ops)[op]) *has = PETSC_TRUE; 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M) 118 { 119 Mat_Transpose *Aa = (Mat_Transpose*)A->data; 120 121 PetscFunctionBegin; 122 *M = Aa->A; 123 PetscFunctionReturn(0); 124 } 125 126 /*@ 127 MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT 128 129 Logically collective on Mat 130 131 Input Parameter: 132 . A - the MATTRANSPOSE matrix 133 134 Output Parameter: 135 . M - the matrix object stored inside A 136 137 Level: intermediate 138 139 .seealso: MatCreateTranspose() 140 141 @*/ 142 PetscErrorCode MatTransposeGetMat(Mat A,Mat *M) 143 { 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 148 PetscValidType(A,1); 149 PetscValidPointer(M,2); 150 ierr = PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));CHKERRQ(ierr); 151 PetscFunctionReturn(0); 152 } 153 154 /*@ 155 MatCreateTranspose - Creates a new matrix object that behaves like A' 156 157 Collective on Mat 158 159 Input Parameter: 160 . A - the (possibly rectangular) matrix 161 162 Output Parameter: 163 . N - the matrix that represents A' 164 165 Level: intermediate 166 167 Notes: 168 The transpose A' is NOT actually formed! Rather the new matrix 169 object performs the matrix-vector product by using the MatMultTranspose() on 170 the original matrix 171 172 .seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate() 173 174 @*/ 175 PetscErrorCode MatCreateTranspose(Mat A,Mat *N) 176 { 177 PetscErrorCode ierr; 178 PetscInt m,n; 179 Mat_Transpose *Na; 180 181 PetscFunctionBegin; 182 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 183 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr); 184 ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 185 ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr); 186 ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr); 187 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr); 188 189 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 190 (*N)->data = (void*) Na; 191 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 192 Na->A = A; 193 194 (*N)->ops->destroy = MatDestroy_Transpose; 195 (*N)->ops->mult = MatMult_Transpose; 196 (*N)->ops->multadd = MatMultAdd_Transpose; 197 (*N)->ops->multtranspose = MatMultTranspose_Transpose; 198 (*N)->ops->multtransposeadd = MatMultTransposeAdd_Transpose; 199 (*N)->ops->duplicate = MatDuplicate_Transpose; 200 (*N)->ops->getvecs = MatCreateVecs_Transpose; 201 (*N)->ops->axpy = MatAXPY_Transpose; 202 (*N)->ops->hasoperation = MatHasOperation_Transpose; 203 (*N)->assembled = PETSC_TRUE; 204 205 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);CHKERRQ(ierr); 206 ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 207 ierr = MatSetUp(*N);CHKERRQ(ierr); 208 PetscFunctionReturn(0); 209 } 210 211