Actual source code: normmh.c

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

      4: typedef struct {
      5:   Mat         A;
      6:   Vec         w,left,right,leftwork,rightwork;
      7:   PetscScalar scale;
      8: } Mat_Normal;

     10: PetscErrorCode MatScaleHermitian_Normal(Mat inA,PetscScalar scale)
     11: {
     12:   Mat_Normal *a = (Mat_Normal*)inA->data;

     15:   a->scale *= scale;
     16:   return(0);
     17: }

     19: PetscErrorCode MatDiagonalScaleHermitian_Normal(Mat inA,Vec left,Vec right)
     20: {
     21:   Mat_Normal     *a = (Mat_Normal*)inA->data;

     25:   if (left) {
     26:     if (!a->left) {
     27:       VecDuplicate(left,&a->left);
     28:       VecCopy(left,a->left);
     29:     } else {
     30:       VecPointwiseMult(a->left,left,a->left);
     31:     }
     32:   }
     33:   if (right) {
     34:     if (!a->right) {
     35:       VecDuplicate(right,&a->right);
     36:       VecCopy(right,a->right);
     37:     } else {
     38:       VecPointwiseMult(a->right,right,a->right);
     39:     }
     40:   }
     41:   return(0);
     42: }

     44: PetscErrorCode MatMultHermitian_Normal(Mat N,Vec x,Vec y)
     45: {
     46:   Mat_Normal     *Na = (Mat_Normal*)N->data;
     48:   Vec            in;

     51:   in = x;
     52:   if (Na->right) {
     53:     if (!Na->rightwork) {
     54:       VecDuplicate(Na->right,&Na->rightwork);
     55:     }
     56:     VecPointwiseMult(Na->rightwork,Na->right,in);
     57:     in   = Na->rightwork;
     58:   }
     59:   MatMult(Na->A,in,Na->w);
     60:   MatMultHermitianTranspose(Na->A,Na->w,y);
     61:   if (Na->left) {
     62:     VecPointwiseMult(y,Na->left,y);
     63:   }
     64:   VecScale(y,Na->scale);
     65:   return(0);
     66: }

     68: PetscErrorCode MatMultHermitianAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
     69: {
     70:   Mat_Normal     *Na = (Mat_Normal*)N->data;
     72:   Vec            in;

     75:   in = v1;
     76:   if (Na->right) {
     77:     if (!Na->rightwork) {
     78:       VecDuplicate(Na->right,&Na->rightwork);
     79:     }
     80:     VecPointwiseMult(Na->rightwork,Na->right,in);
     81:     in   = Na->rightwork;
     82:   }
     83:   MatMult(Na->A,in,Na->w);
     84:   VecScale(Na->w,Na->scale);
     85:   if (Na->left) {
     86:     MatMultHermitianTranspose(Na->A,Na->w,v3);
     87:     VecPointwiseMult(v3,Na->left,v3);
     88:     VecAXPY(v3,1.0,v2);
     89:   } else {
     90:     MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);
     91:   }
     92:   return(0);
     93: }

     95: PetscErrorCode MatMultHermitianTranspose_Normal(Mat N,Vec x,Vec y)
     96: {
     97:   Mat_Normal     *Na = (Mat_Normal*)N->data;
     99:   Vec            in;

    102:   in = x;
    103:   if (Na->left) {
    104:     if (!Na->leftwork) {
    105:       VecDuplicate(Na->left,&Na->leftwork);
    106:     }
    107:     VecPointwiseMult(Na->leftwork,Na->left,in);
    108:     in   = Na->leftwork;
    109:   }
    110:   MatMult(Na->A,in,Na->w);
    111:   MatMultHermitianTranspose(Na->A,Na->w,y);
    112:   if (Na->right) {
    113:     VecPointwiseMult(y,Na->right,y);
    114:   }
    115:   VecScale(y,Na->scale);
    116:   return(0);
    117: }

    119: PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
    120: {
    121:   Mat_Normal     *Na = (Mat_Normal*)N->data;
    123:   Vec            in;

    126:   in = v1;
    127:   if (Na->left) {
    128:     if (!Na->leftwork) {
    129:       VecDuplicate(Na->left,&Na->leftwork);
    130:     }
    131:     VecPointwiseMult(Na->leftwork,Na->left,in);
    132:     in   = Na->leftwork;
    133:   }
    134:   MatMult(Na->A,in,Na->w);
    135:   VecScale(Na->w,Na->scale);
    136:   if (Na->right) {
    137:     MatMultHermitianTranspose(Na->A,Na->w,v3);
    138:     VecPointwiseMult(v3,Na->right,v3);
    139:     VecAXPY(v3,1.0,v2);
    140:   } else {
    141:     MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);
    142:   }
    143:   return(0);
    144: }

    146: PetscErrorCode MatDestroyHermitian_Normal(Mat N)
    147: {
    148:   Mat_Normal     *Na = (Mat_Normal*)N->data;

    152:   MatDestroy(&Na->A);
    153:   VecDestroy(&Na->w);
    154:   VecDestroy(&Na->left);
    155:   VecDestroy(&Na->right);
    156:   VecDestroy(&Na->leftwork);
    157:   VecDestroy(&Na->rightwork);
    158:   PetscFree(N->data);
    159:   return(0);
    160: }

    162: /*
    163:       Slow, nonscalable version
    164: */
    165: PetscErrorCode MatGetDiagonalHermitian_Normal(Mat N,Vec v)
    166: {
    167:   Mat_Normal        *Na = (Mat_Normal*)N->data;
    168:   Mat               A   = Na->A;
    169:   PetscErrorCode    ierr;
    170:   PetscInt          i,j,rstart,rend,nnz;
    171:   const PetscInt    *cols;
    172:   PetscScalar       *diag,*work,*values;
    173:   const PetscScalar *mvalues;

    176:   PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);
    177:   PetscArrayzero(work,A->cmap->N);
    178:   MatGetOwnershipRange(A,&rstart,&rend);
    179:   for (i=rstart; i<rend; i++) {
    180:     MatGetRow(A,i,&nnz,&cols,&mvalues);
    181:     for (j=0; j<nnz; j++) {
    182:       work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]);
    183:     }
    184:     MatRestoreRow(A,i,&nnz,&cols,&mvalues);
    185:   }
    186:   MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
    187:   rstart = N->cmap->rstart;
    188:   rend   = N->cmap->rend;
    189:   VecGetArray(v,&values);
    190:   PetscArraycpy(values,diag+rstart,rend-rstart);
    191:   VecRestoreArray(v,&values);
    192:   PetscFree2(diag,work);
    193:   VecScale(v,Na->scale);
    194:   return(0);
    195: }

    197: /*@
    198:       MatCreateNormalHermitian - Creates a new matrix object that behaves like (A*)'*A.

    200:    Collective on Mat

    202:    Input Parameter:
    203: .   A  - the (possibly rectangular complex) matrix

    205:    Output Parameter:
    206: .   N - the matrix that represents (A*)'*A

    208:    Level: intermediate

    210:    Notes:
    211:     The product (A*)'*A is NOT actually formed! Rather the new matrix
    212:           object performs the matrix-vector product by first multiplying by
    213:           A and then (A*)'
    214: @*/
    215: PetscErrorCode  MatCreateNormalHermitian(Mat A,Mat *N)
    216: {
    218:   PetscInt       m,n;
    219:   Mat_Normal     *Na;

    222:   MatGetLocalSize(A,&m,&n);
    223:   MatCreate(PetscObjectComm((PetscObject)A),N);
    224:   MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
    225:   PetscObjectChangeTypeName((PetscObject)*N,MATNORMALHERMITIAN);

    227:   PetscNewLog(*N,&Na);
    228:   (*N)->data = (void*) Na;
    229:   PetscObjectReference((PetscObject)A);
    230:   Na->A      = A;
    231:   Na->scale  = 1.0;

    233:   VecCreateMPI(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,&Na->w);

    235:   (*N)->ops->destroy          = MatDestroyHermitian_Normal;
    236:   (*N)->ops->mult             = MatMultHermitian_Normal;
    237:   (*N)->ops->multtranspose    = MatMultHermitianTranspose_Normal;
    238:   (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal;
    239:   (*N)->ops->multadd          = MatMultHermitianAdd_Normal;
    240:   (*N)->ops->getdiagonal      = MatGetDiagonalHermitian_Normal;
    241:   (*N)->ops->scale            = MatScaleHermitian_Normal;
    242:   (*N)->ops->diagonalscale    = MatDiagonalScaleHermitian_Normal;
    243:   (*N)->assembled             = PETSC_TRUE;
    244:   (*N)->cmap->N               = A->cmap->N;
    245:   (*N)->rmap->N               = A->cmap->N;
    246:   (*N)->cmap->n               = A->cmap->n;
    247:   (*N)->rmap->n               = A->cmap->n;

    249:   MatSetOption(*N,MAT_HERMITIAN,PETSC_TRUE);
    250:   return(0);
    251: }