Actual source code: htransm.c

      2: #include <petsc/private/matimpl.h>

      4: typedef struct {
      5:   Mat A;
      6: } Mat_HT;

      8: PetscErrorCode MatMult_HT(Mat N,Vec x,Vec y)
      9: {
     10:   Mat_HT         *Na = (Mat_HT*)N->data;

     14:   MatMultHermitianTranspose(Na->A,x,y);
     15:   return(0);
     16: }

     18: PetscErrorCode MatMultAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)
     19: {
     20:   Mat_HT         *Na = (Mat_HT*)N->data;

     24:   MatMultHermitianTransposeAdd(Na->A,v1,v2,v3);
     25:   return(0);
     26: }

     28: PetscErrorCode MatMultHermitianTranspose_HT(Mat N,Vec x,Vec y)
     29: {
     30:   Mat_HT         *Na = (Mat_HT*)N->data;

     34:   MatMult(Na->A,x,y);
     35:   return(0);
     36: }

     38: PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)
     39: {
     40:   Mat_HT         *Na = (Mat_HT*)N->data;

     44:   MatMultAdd(Na->A,v1,v2,v3);
     45:   return(0);
     46: }

     48: PetscErrorCode MatDestroy_HT(Mat N)
     49: {
     50:   Mat_HT         *Na = (Mat_HT*)N->data;

     54:   MatDestroy(&Na->A);
     55:   PetscObjectComposeFunction((PetscObject)N,"MatHermitianTransposeGetMat_C",NULL);
     56: #if !defined(PETSC_USE_COMPLEX)
     57:   PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
     58:   PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_anytype_C",NULL);
     59: #endif
     60:   PetscFree(N->data);
     61:   return(0);
     62: }

     64: PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat* m)
     65: {
     66:   Mat_HT         *Na = (Mat_HT*)N->data;

     70:   if (op == MAT_COPY_VALUES) {
     71:     MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,m);
     72:   } else if (op == MAT_DO_NOT_COPY_VALUES) {
     73:     MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);
     74:     MatHermitianTranspose(*m,MAT_INPLACE_MATRIX,m);
     75:   } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
     76:   return(0);
     77: }

     79: PetscErrorCode MatCreateVecs_HT(Mat N,Vec *r, Vec *l)
     80: {
     81:   Mat_HT         *Na = (Mat_HT*)N->data;

     85:   MatCreateVecs(Na->A,l,r);
     86:   return(0);
     87: }

     89: PetscErrorCode MatAXPY_HT(Mat Y,PetscScalar a,Mat X,MatStructure str)
     90: {
     91:   Mat_HT         *Ya = (Mat_HT*)Y->data;
     92:   Mat_HT         *Xa = (Mat_HT*)X->data;
     93:   Mat              M = Ya->A;
     94:   Mat              N = Xa->A;

     98:   MatAXPY(M,a,N,str);
     99:   return(0);
    100: }

    102: PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N,Mat *M)
    103: {
    104:   Mat_HT         *Na = (Mat_HT*)N->data;

    107:   *M = Na->A;
    108:   return(0);
    109: }

    111: /*@
    112:       MatHermitianTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT

    114:    Logically collective on Mat

    116:    Input Parameter:
    117: .   A  - the MATTRANSPOSE matrix

    119:    Output Parameter:
    120: .   M - the matrix object stored inside A

    122:    Level: intermediate

    124: .seealso: MatCreateHermitianTranspose()

    126: @*/
    127: PetscErrorCode MatHermitianTransposeGetMat(Mat A,Mat *M)
    128: {

    135:   PetscUseMethod(A,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(A,M));
    136:   return(0);
    137: }

    139: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat);

    141: PetscErrorCode MatGetDiagonal_HT(Mat A,Vec v)
    142: {
    143:   Mat_HT         *Na = (Mat_HT*)A->data;

    147:   MatGetDiagonal(Na->A,v);
    148:   VecConjugate(v);
    149:   return(0);
    150: }

    152: PetscErrorCode MatConvert_HT(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
    153: {
    154:   Mat_HT         *Na = (Mat_HT*)A->data;
    155:   Mat            B;

    159:   MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,&B);
    160:   if (reuse != MAT_INPLACE_MATRIX) {
    161:     MatConvert(B,newtype,reuse,newmat);
    162:     MatDestroy(&B);
    163:   } else {
    164:     MatConvert(B,newtype,MAT_INPLACE_MATRIX,&B);
    165:     MatHeaderReplace(A,&B);
    166:   }
    167:   return(0);
    168: }

    170: /*@
    171:       MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'*

    173:    Collective on Mat

    175:    Input Parameter:
    176: .   A  - the (possibly rectangular) matrix

    178:    Output Parameter:
    179: .   N - the matrix that represents A'*

    181:    Level: intermediate

    183:    Notes:
    184:     The hermitian transpose A' is NOT actually formed! Rather the new matrix
    185:           object performs the matrix-vector product by using the MatMultHermitianTranspose() on
    186:           the original matrix

    188: .seealso: MatCreateNormal(), MatMult(), MatMultHermitianTranspose(), MatCreate()

    190: @*/
    191: PetscErrorCode  MatCreateHermitianTranspose(Mat A,Mat *N)
    192: {
    194:   PetscInt       m,n;
    195:   Mat_HT         *Na;
    196:   VecType        vtype;

    199:   MatGetLocalSize(A,&m,&n);
    200:   MatCreate(PetscObjectComm((PetscObject)A),N);
    201:   MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
    202:   PetscLayoutSetUp((*N)->rmap);
    203:   PetscLayoutSetUp((*N)->cmap);
    204:   PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);

    206:   PetscNewLog(*N,&Na);
    207:   (*N)->data = (void*) Na;
    208:   PetscObjectReference((PetscObject)A);
    209:   Na->A      = A;

    211:   (*N)->ops->destroy                   = MatDestroy_HT;
    212:   (*N)->ops->mult                      = MatMult_HT;
    213:   (*N)->ops->multadd                   = MatMultAdd_HT;
    214:   (*N)->ops->multhermitiantranspose    = MatMultHermitianTranspose_HT;
    215:   (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT;
    216:   (*N)->ops->duplicate                 = MatDuplicate_HT;
    217:   (*N)->ops->getvecs                   = MatCreateVecs_HT;
    218:   (*N)->ops->axpy                      = MatAXPY_HT;
    219: #if !defined(PETSC_USE_COMPLEX)
    220:   (*N)->ops->productsetfromoptions     = MatProductSetFromOptions_Transpose;
    221: #endif
    222:   (*N)->ops->getdiagonal               = MatGetDiagonal_HT;
    223:   (*N)->ops->convert                   = MatConvert_HT;
    224:   (*N)->assembled                      = PETSC_TRUE;

    226:   PetscObjectComposeFunction((PetscObject)(*N),"MatHermitianTransposeGetMat_C",MatHermitianTransposeGetMat_HT);
    227: #if !defined(PETSC_USE_COMPLEX)
    228:   PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatHermitianTransposeGetMat_HT);
    229:   PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_anytype_C",MatProductSetFromOptions_Transpose);
    230: #endif
    231:   MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
    232:   MatGetVecType(A,&vtype);
    233:   MatSetVecType(*N,vtype);
    234: #if defined(PETSC_HAVE_DEVICE)
    235:   MatBindToCPU(*N,A->boundtocpu);
    236: #endif
    237:   MatSetUp(*N);
    238:   return(0);
    239: }