Actual source code: axpy.c

      1: #define PETSCMAT_DLL

     3:  #include src/mat/matimpl.h

      7: /*@
      8:    MatAXPY - Computes Y = a*X + Y.

     10:    Collective on Mat

     12:    Input Parameters:
     13: +  a - the scalar multiplier
     14: .  X - the first matrix
     15: .  Y - the second matrix
     16: -  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN

     18:    Contributed by: Matthew Knepley

     20:    Notes:
     21:      Will only be efficient if one has the SAME_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN

     23:    Level: intermediate

     25: .keywords: matrix, add

     27: .seealso: MatAYPX()
     28:  @*/
     29: PetscErrorCode PETSCMAT_DLLEXPORT MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str)
     30: {
     32:   PetscInt       m1,m2,n1,n2;

     38:   MatGetSize(X,&m1,&n1);
     39:   MatGetSize(Y,&m2,&n2);
     40:   if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrix add: %D %D %D %D",m1,m2,n1,n2);

     42:   if (Y->ops->axpy) {
     43:     (*Y->ops->axpy)(Y,a,X,str);
     44:   } else {
     45:     MatAXPY_Basic(Y,a,X,str);
     46:   }
     47:   return(0);
     48: }

     53: PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str)
     54: {
     55:   PetscInt          i,start,end,j,ncols,m,n;
     56:   PetscErrorCode    ierr;
     57:   const PetscInt    *row;
     58:   PetscScalar       *val;
     59:   const PetscScalar *vals;

     62:   MatGetSize(X,&m,&n);
     63:   MatGetOwnershipRange(X,&start,&end);
     64:   if (a == 1.0) {
     65:     for (i = start; i < end; i++) {
     66:       MatGetRow(X,i,&ncols,&row,&vals);
     67:       MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);
     68:       MatRestoreRow(X,i,&ncols,&row,&vals);
     69:     }
     70:   } else {
     71:     PetscMalloc((n+1)*sizeof(PetscScalar),&val);
     72:     for (i=start; i<end; i++) {
     73:       MatGetRow(X,i,&ncols,&row,&vals);
     74:       for (j=0; j<ncols; j++) {
     75:         val[j] = a*vals[j];
     76:       }
     77:       MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);
     78:       MatRestoreRow(X,i,&ncols,&row,&vals);
     79:     }
     80:     PetscFree(val);
     81:   }
     82:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
     83:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
     84:   return(0);
     85: }

     89: /*@
     90:    MatShift - Computes Y =  Y + a I, where a is a PetscScalar and I is the identity matrix.

     92:    Collective on Mat

     94:    Input Parameters:
     95: +  Y - the matrices
     96: -  a - the PetscScalar

     98:    Level: intermediate

    100: .keywords: matrix, add, shift

    102: .seealso: MatDiagonalSet()
    103:  @*/
    104: PetscErrorCode PETSCMAT_DLLEXPORT MatShift(Mat Y,PetscScalar a)
    105: {
    107:   PetscInt       i,start,end;

    111:   if (!Y->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    112:   if (Y->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    113:   MatPreallocated(Y);

    115:   if (Y->ops->shift) {
    116:     (*Y->ops->shift)(Y,a);
    117:   } else {
    118:     PetscScalar alpha = a;
    119:     MatGetOwnershipRange(Y,&start,&end);
    120:     for (i=start; i<end; i++) {
    121:       MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);
    122:     }
    123:     MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
    124:     MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
    125:   }
    126:   return(0);
    127: }

    131: /*@
    132:    MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix
    133:    that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is
    134:    INSERT_VALUES.

    136:    Input Parameters:
    137: +  Y - the input matrix
    138: .  D - the diagonal matrix, represented as a vector
    139: -  i - INSERT_VALUES or ADD_VALUES

    141:    Collective on Mat and Vec

    143:    Level: intermediate

    145: .keywords: matrix, add, shift, diagonal

    147: .seealso: MatShift()
    148: @*/
    149: PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalSet(Mat Y,Vec D,InsertMode is)
    150: {
    152:   PetscInt       i,start,end;

    157:   if (Y->ops->diagonalset) {
    158:     (*Y->ops->diagonalset)(Y,D,is);
    159:   } else {
    160:     PetscInt    vstart,vend;
    161:     PetscScalar *v;
    162:     VecGetOwnershipRange(D,&vstart,&vend);
    163:     MatGetOwnershipRange(Y,&start,&end);
    164:     if (vstart != start || vend != end) {
    165:       SETERRQ4(PETSC_ERR_ARG_SIZ,"Vector ownership range not compatible with matrix: %D %D vec %D %D mat",vstart,vend,start,end);
    166:     }
    167:     VecGetArray(D,&v);
    168:     for (i=start; i<end; i++) {
    169:       MatSetValues(Y,1,&i,1,&i,v+i-start,is);
    170:     }
    171:     VecRestoreArray(D,&v);
    172:     MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
    173:     MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
    174:   }
    175:   return(0);
    176: }

    180: /*@
    181:    MatAYPX - Computes Y = X + a*Y.

    183:    Collective on Mat

    185:    Input Parameters:
    186: +  X,Y - the matrices
    187: -  a - the PetscScalar multiplier

    189:    Contributed by: Matthew Knepley

    191:    Notes:
    192:    This routine currently uses the MatAXPY() implementation.

    194:    This is slow, if you need it fast send email to petsc-maint@mcs.anl.gov

    196:    Level: intermediate

    198: .keywords: matrix, add

    200: .seealso: MatAXPY()
    201:  @*/
    202: PetscErrorCode PETSCMAT_DLLEXPORT MatAYPX(Mat Y,PetscScalar a,Mat X)
    203: {
    204:   PetscScalar    one = 1.0;
    206:   PetscInt       mX,mY,nX,nY;

    212:   MatGetSize(X,&mX,&nX);
    213:   MatGetSize(X,&mY,&nY);
    214:   if (mX != mY || nX != nY) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrices: %D %D first %D %D second",mX,mY,nX,nY);

    216:   MatScale(Y,a);
    217:   MatAXPY(Y,one,X,DIFFERENT_NONZERO_PATTERN);
    218:   return(0);
    219: }

    223: /*@
    224:     MatComputeExplicitOperator - Computes the explicit matrix

    226:     Collective on Mat

    228:     Input Parameter:
    229: .   inmat - the matrix

    231:     Output Parameter:
    232: .   mat - the explict preconditioned operator

    234:     Notes:
    235:     This computation is done by applying the operators to columns of the
    236:     identity matrix.

    238:     Currently, this routine uses a dense matrix format when 1 processor
    239:     is used and a sparse format otherwise.  This routine is costly in general,
    240:     and is recommended for use only with relatively small systems.

    242:     Level: advanced
    243:
    244: .keywords: Mat, compute, explicit, operator

    246: @*/
    247: PetscErrorCode PETSCMAT_DLLEXPORT MatComputeExplicitOperator(Mat inmat,Mat *mat)
    248: {
    249:   Vec            in,out;
    251:   PetscInt       i,M,m,*rows,start,end;
    252:   MPI_Comm       comm;
    253:   PetscScalar    *array,zero = 0.0,one = 1.0;
    254:   PetscMPIInt    size;

    260:   comm = inmat->comm;
    261:   MPI_Comm_size(comm,&size);

    263:   MatGetLocalSize(inmat,&m,0);
    264:   MatGetSize(inmat,&M,0);
    265:   VecCreateMPI(comm,m,M,&in);
    266:   VecDuplicate(in,&out);
    267:   VecGetOwnershipRange(in,&start,&end);
    268:   PetscMalloc((m+1)*sizeof(PetscInt),&rows);
    269:   for (i=0; i<m; i++) {rows[i] = start + i;}

    271:   MatCreate(comm,mat);
    272:   MatSetSizes(*mat,m,m,M,M);
    273:   if (size == 1) {
    274:     MatSetType(*mat,MATSEQDENSE);
    275:     MatSeqDenseSetPreallocation(*mat,PETSC_NULL);
    276:   } else {
    277:     MatSetType(*mat,MATMPIAIJ);
    278:     MatMPIAIJSetPreallocation(*mat,0,PETSC_NULL,0,PETSC_NULL);
    279:   }

    281:   for (i=0; i<M; i++) {

    283:     VecSet(in,zero);
    284:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
    285:     VecAssemblyBegin(in);
    286:     VecAssemblyEnd(in);

    288:     MatMult(inmat,in,out);
    289:
    290:     VecGetArray(out,&array);
    291:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
    292:     VecRestoreArray(out,&array);

    294:   }
    295:   PetscFree(rows);
    296:   VecDestroy(out);
    297:   VecDestroy(in);
    298:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
    299:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
    300:   return(0);
    301: }

    303: /* Get the map xtoy which is used by MatAXPY() in the case of SUBSET_NONZERO_PATTERN */
    306: PetscErrorCode MatAXPYGetxtoy_Private(PetscInt m,PetscInt *xi,PetscInt *xj,PetscInt *xgarray, PetscInt *yi,PetscInt *yj,PetscInt *ygarray, PetscInt **xtoy)
    307: {
    309:   PetscInt       row,i,nz,xcol,ycol,jx,jy,*x2y;

    312:   PetscMalloc(xi[m]*sizeof(PetscInt),&x2y);
    313:   i = 0;
    314:   for (row=0; row<m; row++){
    315:     nz = xi[1] - xi[0];
    316:     jy = 0;
    317:     for (jx=0; jx<nz; jx++,jy++){
    318:       if (xgarray && ygarray){
    319:         xcol = xgarray[xj[*xi + jx]];
    320:         ycol = ygarray[yj[*yi + jy]];
    321:       } else {
    322:         xcol = xj[*xi + jx];
    323:         ycol = yj[*yi + jy];  /* col index for y */
    324:       }
    325:       while ( ycol < xcol ) {
    326:         jy++;
    327:         if (ygarray){
    328:           ycol = ygarray[yj[*yi + jy]];
    329:         } else {
    330:           ycol = yj[*yi + jy];
    331:         }
    332:       }
    333:       if (xcol != ycol) SETERRQ2(PETSC_ERR_ARG_WRONG,"X matrix entry (%D,%D) is not in Y matrix",row,ycol);
    334:       x2y[i++] = *yi + jy;
    335:     }
    336:     xi++; yi++;
    337:   }
    338:   *xtoy = x2y;
    339:   return(0);
    340: }