Actual source code: dense.c

      1: #define PETSCMAT_DLL

      3: /*
      4:      Defines the basic matrix operations for sequential dense.
      5: */

     7:  #include src/mat/impls/dense/seq/dense.h
     8:  #include petscblaslapack.h

     12: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
     13: {
     14:   Mat_SeqDense   *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
     15:   PetscScalar    oalpha = alpha;
     16:   PetscInt       j;
     17:   PetscBLASInt   N = (PetscBLASInt)X->m*X->n,m=(PetscBLASInt)X->m,ldax = x->lda,lday=y->lda,one = 1;

     21:   if (X->m != Y->m || X->n != Y->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(X) != size(Y)");
     22:   if (ldax>m || lday>m) {
     23:     for (j=0; j<X->n; j++) {
     24:       BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one);
     25:     }
     26:   } else {
     27:     BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one);
     28:   }
     29:   PetscLogFlops(2*N-1);
     30:   return(0);
     31: }

     35: PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
     36: {
     37:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
     38:   PetscInt     i,N = A->m*A->n,count = 0;
     39:   PetscScalar  *v = a->v;

     42:   for (i=0; i<N; i++) {if (*v != 0.0) count++; v++;}

     44:   info->rows_global       = (double)A->m;
     45:   info->columns_global    = (double)A->n;
     46:   info->rows_local        = (double)A->m;
     47:   info->columns_local     = (double)A->n;
     48:   info->block_size        = 1.0;
     49:   info->nz_allocated      = (double)N;
     50:   info->nz_used           = (double)count;
     51:   info->nz_unneeded       = (double)(N-count);
     52:   info->assemblies        = (double)A->num_ass;
     53:   info->mallocs           = 0;
     54:   info->memory            = A->mem;
     55:   info->fill_ratio_given  = 0;
     56:   info->fill_ratio_needed = 0;
     57:   info->factor_mallocs    = 0;

     59:   return(0);
     60: }

     64: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
     65: {
     66:   Mat_SeqDense *a = (Mat_SeqDense*)A->data;
     67:   PetscScalar  oalpha = alpha;
     68:   PetscBLASInt one = 1,lda = a->lda,j,nz;

     72:   if (lda>A->m) {
     73:     nz = (PetscBLASInt)A->m;
     74:     for (j=0; j<A->n; j++) {
     75:       BLASscal_(&nz,&oalpha,a->v+j*lda,&one);
     76:     }
     77:   } else {
     78:     nz = (PetscBLASInt)A->m*A->n;
     79:     BLASscal_(&nz,&oalpha,a->v,&one);
     80:   }
     81:   PetscLogFlops(nz);
     82:   return(0);
     83: }
     84:
     85: /* ---------------------------------------------------------------*/
     86: /* COMMENT: I have chosen to hide row permutation in the pivots,
     87:    rather than put it in the Mat->row slot.*/
     90: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,MatFactorInfo *minfo)
     91: {
     92: #if defined(PETSC_MISSING_LAPACK_GETRF)
     94:   SETERRQ(PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
     95: #else
     96:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
     98:   PetscBLASInt   n = (PetscBLASInt)A->n,m = (PetscBLASInt)A->m,info;

    101:   if (!mat->pivots) {
    102:     PetscMalloc((A->m+1)*sizeof(PetscBLASInt),&mat->pivots);
    103:     PetscLogObjectMemory(A,A->m*sizeof(PetscBLASInt));
    104:   }
    105:   A->factor = FACTOR_LU;
    106:   if (!A->m || !A->n) return(0);
    107:   LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info);
    108:   if (info<0) SETERRQ(PETSC_ERR_LIB,"Bad argument to LU factorization");
    109:   if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
    110:   PetscLogFlops((2*A->n*A->n*A->n)/3);
    111: #endif
    112:   return(0);
    113: }

    117: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
    118: {
    119:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l;
    121:   PetscInt       lda = (PetscInt)mat->lda,j,m;
    122:   Mat            newi;

    125:   MatCreate(A->comm,&newi);
    126:   MatSetSizes(newi,A->m,A->n,A->m,A->n);
    127:   MatSetType(newi,A->type_name);
    128:   MatSeqDenseSetPreallocation(newi,PETSC_NULL);
    129:   if (cpvalues == MAT_COPY_VALUES) {
    130:     l = (Mat_SeqDense*)newi->data;
    131:     if (lda>A->m) {
    132:       m = A->m;
    133:       for (j=0; j<A->n; j++) {
    134:         PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
    135:       }
    136:     } else {
    137:       PetscMemcpy(l->v,mat->v,A->m*A->n*sizeof(PetscScalar));
    138:     }
    139:   }
    140:   newi->assembled = PETSC_TRUE;
    141:   *newmat = newi;
    142:   return(0);
    143: }

    147: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,MatFactorInfo *info,Mat *fact)
    148: {

    152:   MatDuplicate_SeqDense(A,MAT_DO_NOT_COPY_VALUES,fact);
    153:   return(0);
    154: }

    158: PetscErrorCode MatLUFactorNumeric_SeqDense(Mat A,MatFactorInfo *info_dummy,Mat *fact)
    159: {
    160:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data,*l = (Mat_SeqDense*)(*fact)->data;
    162:   PetscInt       lda1=mat->lda,lda2=l->lda, m=A->m,n=A->n, j;
    163:   MatFactorInfo  info;

    166:   /* copy the numerical values */
    167:   if (lda1>m || lda2>m ) {
    168:     for (j=0; j<n; j++) {
    169:       PetscMemcpy(l->v+j*lda2,mat->v+j*lda1,m*sizeof(PetscScalar));
    170:     }
    171:   } else {
    172:     PetscMemcpy(l->v,mat->v,A->m*A->n*sizeof(PetscScalar));
    173:   }
    174:   (*fact)->factor = 0;
    175:   MatLUFactor(*fact,0,0,&info);
    176:   return(0);
    177: }

    181: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,MatFactorInfo *info,Mat *fact)
    182: {

    186:   MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,fact);
    187:   return(0);
    188: }

    192: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,MatFactorInfo *factinfo)
    193: {
    194: #if defined(PETSC_MISSING_LAPACK_POTRF)
    196:   SETERRQ(PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
    197: #else
    198:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    200:   PetscBLASInt   n = (PetscBLASInt)A->n,info;
    201:
    203:   if (mat->pivots) {
    204:     PetscFree(mat->pivots);
    205:     PetscLogObjectMemory(A,-A->m*sizeof(PetscInt));
    206:     mat->pivots = 0;
    207:   }
    208:   if (!A->m || !A->n) return(0);
    209:   LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info);
    210:   if (info) SETERRQ1(PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
    211:   A->factor = FACTOR_CHOLESKY;
    212:   PetscLogFlops((A->n*A->n*A->n)/3);
    213: #endif
    214:   return(0);
    215: }

    219: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat A,MatFactorInfo *info_dummy,Mat *fact)
    220: {
    222:   MatFactorInfo  info;

    225:   info.fill = 1.0;
    226:   MatCholeskyFactor_SeqDense(*fact,0,&info);
    227:   return(0);
    228: }

    232: PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
    233: {
    234:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    236:   PetscBLASInt   m = (PetscBLASInt)A->m, one = 1,info;
    237:   PetscScalar    *x,*y;
    238:
    240:   if (!A->m || !A->n) return(0);
    241:   VecGetArray(xx,&x);
    242:   VecGetArray(yy,&y);
    243:   PetscMemcpy(y,x,A->m*sizeof(PetscScalar));
    244:   if (A->factor == FACTOR_LU) {
    245: #if defined(PETSC_MISSING_LAPACK_GETRS)
    246:     SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
    247: #else
    248:     LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
    249:     if (info) SETERRQ(PETSC_ERR_LIB,"GETRS - Bad solve");
    250: #endif
    251:   } else if (A->factor == FACTOR_CHOLESKY){
    252: #if defined(PETSC_MISSING_LAPACK_POTRS)
    253:     SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
    254: #else
    255:     LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
    256:     if (info) SETERRQ(PETSC_ERR_LIB,"POTRS Bad solve");
    257: #endif
    258:   }
    259:   else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
    260:   VecRestoreArray(xx,&x);
    261:   VecRestoreArray(yy,&y);
    262:   PetscLogFlops(2*A->n*A->n - A->n);
    263:   return(0);
    264: }

    268: PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
    269: {
    270:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    272:   PetscBLASInt   m = (PetscBLASInt) A->m,one = 1,info;
    273:   PetscScalar    *x,*y;
    274:
    276:   if (!A->m || !A->n) return(0);
    277:   VecGetArray(xx,&x);
    278:   VecGetArray(yy,&y);
    279:   PetscMemcpy(y,x,A->m*sizeof(PetscScalar));
    280:   /* assume if pivots exist then use LU; else Cholesky */
    281:   if (mat->pivots) {
    282: #if defined(PETSC_MISSING_LAPACK_GETRS)
    283:     SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
    284: #else
    285:     LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
    286:     if (info) SETERRQ(PETSC_ERR_LIB,"POTRS - Bad solve");
    287: #endif
    288:   } else {
    289: #if defined(PETSC_MISSING_LAPACK_POTRS)
    290:     SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
    291: #else
    292:     LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
    293:     if (info) SETERRQ(PETSC_ERR_LIB,"POTRS - Bad solve");
    294: #endif
    295:   }
    296:   VecRestoreArray(xx,&x);
    297:   VecRestoreArray(yy,&y);
    298:   PetscLogFlops(2*A->n*A->n - A->n);
    299:   return(0);
    300: }

    304: PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
    305: {
    306:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    308:   PetscBLASInt   m = (PetscBLASInt)A->m,one = 1,info;
    309:   PetscScalar    *x,*y,sone = 1.0;
    310:   Vec            tmp = 0;
    311:
    313:   VecGetArray(xx,&x);
    314:   VecGetArray(yy,&y);
    315:   if (!A->m || !A->n) return(0);
    316:   if (yy == zz) {
    317:     VecDuplicate(yy,&tmp);
    318:     PetscLogObjectParent(A,tmp);
    319:     VecCopy(yy,tmp);
    320:   }
    321:   PetscMemcpy(y,x,A->m*sizeof(PetscScalar));
    322:   /* assume if pivots exist then use LU; else Cholesky */
    323:   if (mat->pivots) {
    324: #if defined(PETSC_MISSING_LAPACK_GETRS)
    325:     SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
    326: #else
    327:     LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
    328:     if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
    329: #endif
    330:   } else {
    331: #if defined(PETSC_MISSING_LAPACK_POTRS)
    332:     SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
    333: #else
    334:     LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
    335:     if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
    336: #endif
    337:   }
    338:   if (tmp) {VecAXPY(yy,sone,tmp); VecDestroy(tmp);}
    339:   else     {VecAXPY(yy,sone,zz);}
    340:   VecRestoreArray(xx,&x);
    341:   VecRestoreArray(yy,&y);
    342:   PetscLogFlops(2*A->n*A->n);
    343:   return(0);
    344: }

    348: PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
    349: {
    350:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    352:   PetscBLASInt   m = (PetscBLASInt)A->m,one = 1,info;
    353:   PetscScalar    *x,*y,sone = 1.0;
    354:   Vec            tmp;
    355:
    357:   if (!A->m || !A->n) return(0);
    358:   VecGetArray(xx,&x);
    359:   VecGetArray(yy,&y);
    360:   if (yy == zz) {
    361:     VecDuplicate(yy,&tmp);
    362:     PetscLogObjectParent(A,tmp);
    363:     VecCopy(yy,tmp);
    364:   }
    365:   PetscMemcpy(y,x,A->m*sizeof(PetscScalar));
    366:   /* assume if pivots exist then use LU; else Cholesky */
    367:   if (mat->pivots) {
    368: #if defined(PETSC_MISSING_LAPACK_GETRS)
    369:     SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
    370: #else
    371:     LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info);
    372:     if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
    373: #endif
    374:   } else {
    375: #if defined(PETSC_MISSING_LAPACK_POTRS)
    376:     SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
    377: #else
    378:     LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info);
    379:     if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve");
    380: #endif
    381:   }
    382:   if (tmp) {
    383:     VecAXPY(yy,sone,tmp);
    384:     VecDestroy(tmp);
    385:   } else {
    386:     VecAXPY(yy,sone,zz);
    387:   }
    388:   VecRestoreArray(xx,&x);
    389:   VecRestoreArray(yy,&y);
    390:   PetscLogFlops(2*A->n*A->n);
    391:   return(0);
    392: }
    393: /* ------------------------------------------------------------------*/
    396: PetscErrorCode MatRelax_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
    397: {
    398:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    399:   PetscScalar    *x,*b,*v = mat->v,zero = 0.0,xt;
    401:   PetscInt       m = A->m,i;
    402: #if !defined(PETSC_USE_COMPLEX)
    403:   PetscBLASInt   bm = (PetscBLASInt)m, o = 1;
    404: #endif

    407:   if (flag & SOR_ZERO_INITIAL_GUESS) {
    408:     /* this is a hack fix, should have another version without the second BLASdot */
    409:     VecSet(xx,zero);
    410:   }
    411:   VecGetArray(xx,&x);
    412:   VecGetArray(bb,&b);
    413:   its  = its*lits;
    414:   if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
    415:   while (its--) {
    416:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
    417:       for (i=0; i<m; i++) {
    418: #if defined(PETSC_USE_COMPLEX)
    419:         /* cannot use BLAS dot for complex because compiler/linker is
    420:            not happy about returning a double complex */
    421:         PetscInt         _i;
    422:         PetscScalar sum = b[i];
    423:         for (_i=0; _i<m; _i++) {
    424:           sum -= PetscConj(v[i+_i*m])*x[_i];
    425:         }
    426:         xt = sum;
    427: #else
    428:         xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o);
    429: #endif
    430:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
    431:       }
    432:     }
    433:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
    434:       for (i=m-1; i>=0; i--) {
    435: #if defined(PETSC_USE_COMPLEX)
    436:         /* cannot use BLAS dot for complex because compiler/linker is
    437:            not happy about returning a double complex */
    438:         PetscInt         _i;
    439:         PetscScalar sum = b[i];
    440:         for (_i=0; _i<m; _i++) {
    441:           sum -= PetscConj(v[i+_i*m])*x[_i];
    442:         }
    443:         xt = sum;
    444: #else
    445:         xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o);
    446: #endif
    447:         x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
    448:       }
    449:     }
    450:   }
    451:   VecRestoreArray(bb,&b);
    452:   VecRestoreArray(xx,&x);
    453:   return(0);
    454: }

    456: /* -----------------------------------------------------------------*/
    459: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
    460: {
    461:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    462:   PetscScalar    *v = mat->v,*x,*y;
    464:   PetscBLASInt   m = (PetscBLASInt)A->m, n = (PetscBLASInt)A->n,_One=1;
    465:   PetscScalar    _DOne=1.0,_DZero=0.0;

    468:   if (!A->m || !A->n) return(0);
    469:   VecGetArray(xx,&x);
    470:   VecGetArray(yy,&y);
    471:   BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One);
    472:   VecRestoreArray(xx,&x);
    473:   VecRestoreArray(yy,&y);
    474:   PetscLogFlops(2*A->m*A->n - A->n);
    475:   return(0);
    476: }

    480: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
    481: {
    482:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    483:   PetscScalar    *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0;
    485:   PetscBLASInt   m = (PetscBLASInt)A->m, n = (PetscBLASInt)A->n, _One=1;

    488:   if (!A->m || !A->n) return(0);
    489:   VecGetArray(xx,&x);
    490:   VecGetArray(yy,&y);
    491:   BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One);
    492:   VecRestoreArray(xx,&x);
    493:   VecRestoreArray(yy,&y);
    494:   PetscLogFlops(2*A->m*A->n - A->m);
    495:   return(0);
    496: }

    500: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
    501: {
    502:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    503:   PetscScalar    *v = mat->v,*x,*y,_DOne=1.0;
    505:   PetscBLASInt   m = (PetscBLASInt)A->m, n = (PetscBLASInt)A->n, _One=1;

    508:   if (!A->m || !A->n) return(0);
    509:   if (zz != yy) {VecCopy(zz,yy);}
    510:   VecGetArray(xx,&x);
    511:   VecGetArray(yy,&y);
    512:   BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One);
    513:   VecRestoreArray(xx,&x);
    514:   VecRestoreArray(yy,&y);
    515:   PetscLogFlops(2*A->m*A->n);
    516:   return(0);
    517: }

    521: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
    522: {
    523:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    524:   PetscScalar    *v = mat->v,*x,*y;
    526:   PetscBLASInt   m = (PetscBLASInt)A->m, n = (PetscBLASInt)A->n, _One=1;
    527:   PetscScalar    _DOne=1.0;

    530:   if (!A->m || !A->n) return(0);
    531:   if (zz != yy) {VecCopy(zz,yy);}
    532:   VecGetArray(xx,&x);
    533:   VecGetArray(yy,&y);
    534:   BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One);
    535:   VecRestoreArray(xx,&x);
    536:   VecRestoreArray(yy,&y);
    537:   PetscLogFlops(2*A->m*A->n);
    538:   return(0);
    539: }

    541: /* -----------------------------------------------------------------*/
    544: PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
    545: {
    546:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    547:   PetscScalar    *v;
    549:   PetscInt       i;
    550:
    552:   *ncols = A->n;
    553:   if (cols) {
    554:     PetscMalloc((A->n+1)*sizeof(PetscInt),cols);
    555:     for (i=0; i<A->n; i++) (*cols)[i] = i;
    556:   }
    557:   if (vals) {
    558:     PetscMalloc((A->n+1)*sizeof(PetscScalar),vals);
    559:     v    = mat->v + row;
    560:     for (i=0; i<A->n; i++) {(*vals)[i] = *v; v += mat->lda;}
    561:   }
    562:   return(0);
    563: }

    567: PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
    568: {
    571:   if (cols) {PetscFree(*cols);}
    572:   if (vals) {PetscFree(*vals); }
    573:   return(0);
    574: }
    575: /* ----------------------------------------------------------------*/
    578: PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
    579: {
    580:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
    581:   PetscInt     i,j,idx=0;
    582:
    584:   if (!mat->roworiented) {
    585:     if (addv == INSERT_VALUES) {
    586:       for (j=0; j<n; j++) {
    587:         if (indexn[j] < 0) {idx += m; continue;}
    588: #if defined(PETSC_USE_DEBUG)
    589:         if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->n-1);
    590: #endif
    591:         for (i=0; i<m; i++) {
    592:           if (indexm[i] < 0) {idx++; continue;}
    593: #if defined(PETSC_USE_DEBUG)
    594:           if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->m-1);
    595: #endif
    596:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
    597:         }
    598:       }
    599:     } else {
    600:       for (j=0; j<n; j++) {
    601:         if (indexn[j] < 0) {idx += m; continue;}
    602: #if defined(PETSC_USE_DEBUG)
    603:         if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->n-1);
    604: #endif
    605:         for (i=0; i<m; i++) {
    606:           if (indexm[i] < 0) {idx++; continue;}
    607: #if defined(PETSC_USE_DEBUG)
    608:           if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->m-1);
    609: #endif
    610:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
    611:         }
    612:       }
    613:     }
    614:   } else {
    615:     if (addv == INSERT_VALUES) {
    616:       for (i=0; i<m; i++) {
    617:         if (indexm[i] < 0) { idx += n; continue;}
    618: #if defined(PETSC_USE_DEBUG)
    619:         if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->m-1);
    620: #endif
    621:         for (j=0; j<n; j++) {
    622:           if (indexn[j] < 0) { idx++; continue;}
    623: #if defined(PETSC_USE_DEBUG)
    624:           if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->n-1);
    625: #endif
    626:           mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
    627:         }
    628:       }
    629:     } else {
    630:       for (i=0; i<m; i++) {
    631:         if (indexm[i] < 0) { idx += n; continue;}
    632: #if defined(PETSC_USE_DEBUG)
    633:         if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->m-1);
    634: #endif
    635:         for (j=0; j<n; j++) {
    636:           if (indexn[j] < 0) { idx++; continue;}
    637: #if defined(PETSC_USE_DEBUG)
    638:           if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->n-1);
    639: #endif
    640:           mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
    641:         }
    642:       }
    643:     }
    644:   }
    645:   return(0);
    646: }

    650: PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
    651: {
    652:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
    653:   PetscInt     i,j;
    654:   PetscScalar  *vpt = v;

    657:   /* row-oriented output */
    658:   for (i=0; i<m; i++) {
    659:     for (j=0; j<n; j++) {
    660:       *vpt++ = mat->v[indexn[j]*mat->lda + indexm[i]];
    661:     }
    662:   }
    663:   return(0);
    664: }

    666: /* -----------------------------------------------------------------*/

    668:  #include petscsys.h

    672: PetscErrorCode MatLoad_SeqDense(PetscViewer viewer, MatType type,Mat *A)
    673: {
    674:   Mat_SeqDense   *a;
    675:   Mat            B;
    677:   PetscInt       *scols,i,j,nz,header[4];
    678:   int            fd;
    679:   PetscMPIInt    size;
    680:   PetscInt       *rowlengths = 0,M,N,*cols;
    681:   PetscScalar    *vals,*svals,*v,*w;
    682:   MPI_Comm       comm = ((PetscObject)viewer)->comm;

    685:   MPI_Comm_size(comm,&size);
    686:   if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor");
    687:   PetscViewerBinaryGetDescriptor(viewer,&fd);
    688:   PetscBinaryRead(fd,header,4,PETSC_INT);
    689:   if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
    690:   M = header[1]; N = header[2]; nz = header[3];

    692:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
    693:     MatCreate(comm,A);
    694:     MatSetSizes(*A,M,N,M,N);
    695:     MatSetType(*A,type);
    696:     MatSeqDenseSetPreallocation(*A,PETSC_NULL);
    697:     B    = *A;
    698:     a    = (Mat_SeqDense*)B->data;
    699:     v    = a->v;
    700:     /* Allocate some temp space to read in the values and then flip them
    701:        from row major to column major */
    702:     PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);
    703:     /* read in nonzero values */
    704:     PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
    705:     /* now flip the values and store them in the matrix*/
    706:     for (j=0; j<N; j++) {
    707:       for (i=0; i<M; i++) {
    708:         *v++ =w[i*N+j];
    709:       }
    710:     }
    711:     PetscFree(w);
    712:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    713:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    714:   } else {
    715:     /* read row lengths */
    716:     PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);
    717:     PetscBinaryRead(fd,rowlengths,M,PETSC_INT);

    719:     /* create our matrix */
    720:     MatCreate(comm,A);
    721:     MatSetSizes(*A,M,N,M,N);
    722:     MatSetType(*A,type);
    723:     MatSeqDenseSetPreallocation(*A,PETSC_NULL);
    724:     B = *A;
    725:     a = (Mat_SeqDense*)B->data;
    726:     v = a->v;

    728:     /* read column indices and nonzeros */
    729:     PetscMalloc((nz+1)*sizeof(PetscInt),&scols);
    730:     cols = scols;
    731:     PetscBinaryRead(fd,cols,nz,PETSC_INT);
    732:     PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);
    733:     vals = svals;
    734:     PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);

    736:     /* insert into matrix */
    737:     for (i=0; i<M; i++) {
    738:       for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
    739:       svals += rowlengths[i]; scols += rowlengths[i];
    740:     }
    741:     PetscFree(vals);
    742:     PetscFree(cols);
    743:     PetscFree(rowlengths);

    745:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    746:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    747:   }
    748:   return(0);
    749: }

    751:  #include petscsys.h

    755: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
    756: {
    757:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
    758:   PetscErrorCode    ierr;
    759:   PetscInt          i,j;
    760:   const char        *name;
    761:   PetscScalar       *v;
    762:   PetscViewerFormat format;

    765:   PetscObjectGetName((PetscObject)A,&name);
    766:   PetscViewerGetFormat(viewer,&format);
    767:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    768:     return(0);  /* do nothing for now */
    769:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
    770:     PetscViewerASCIIUseTabs(viewer,PETSC_NO);
    771:     for (i=0; i<A->m; i++) {
    772:       v = a->v + i;
    773:       PetscViewerASCIIPrintf(viewer,"row %D:",i);
    774:       for (j=0; j<A->n; j++) {
    775: #if defined(PETSC_USE_COMPLEX)
    776:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
    777:           PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));
    778:         } else if (PetscRealPart(*v)) {
    779:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,PetscRealPart(*v));
    780:         }
    781: #else
    782:         if (*v) {
    783:           PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,*v);
    784:         }
    785: #endif
    786:         v += a->lda;
    787:       }
    788:       PetscViewerASCIIPrintf(viewer,"\n");
    789:     }
    790:     PetscViewerASCIIUseTabs(viewer,PETSC_YES);
    791:   } else {
    792:     PetscViewerASCIIUseTabs(viewer,PETSC_NO);
    793: #if defined(PETSC_USE_COMPLEX)
    794:     PetscTruth allreal = PETSC_TRUE;
    795:     /* determine if matrix has all real values */
    796:     v = a->v;
    797:     for (i=0; i<A->m*A->n; i++) {
    798:         if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;}
    799:     }
    800: #endif
    801:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
    802:       PetscObjectGetName((PetscObject)A,&name);
    803:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->m,A->n);
    804:       PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->m,A->n);
    805:       PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
    806:     }

    808:     for (i=0; i<A->m; i++) {
    809:       v = a->v + i;
    810:       for (j=0; j<A->n; j++) {
    811: #if defined(PETSC_USE_COMPLEX)
    812:         if (allreal) {
    813:           PetscViewerASCIIPrintf(viewer,"%6.4e ",PetscRealPart(*v));
    814:         } else {
    815:           PetscViewerASCIIPrintf(viewer,"%6.4e + %6.4e i ",PetscRealPart(*v),PetscImaginaryPart(*v));
    816:         }
    817: #else
    818:         PetscViewerASCIIPrintf(viewer,"%6.4e ",*v);
    819: #endif
    820:         v += a->lda;
    821:       }
    822:       PetscViewerASCIIPrintf(viewer,"\n");
    823:     }
    824:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
    825:       PetscViewerASCIIPrintf(viewer,"];\n");
    826:     }
    827:     PetscViewerASCIIUseTabs(viewer,PETSC_YES);
    828:   }
    829:   PetscViewerFlush(viewer);
    830:   return(0);
    831: }

    835: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
    836: {
    837:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
    838:   PetscErrorCode    ierr;
    839:   int               fd;
    840:   PetscInt          ict,j,n = A->n,m = A->m,i,*col_lens,nz = m*n;
    841:   PetscScalar       *v,*anonz,*vals;
    842:   PetscViewerFormat format;
    843:
    845:   PetscViewerBinaryGetDescriptor(viewer,&fd);

    847:   PetscViewerGetFormat(viewer,&format);
    848:   if (format == PETSC_VIEWER_BINARY_NATIVE) {
    849:     /* store the matrix as a dense matrix */
    850:     PetscMalloc(4*sizeof(PetscInt),&col_lens);
    851:     col_lens[0] = MAT_FILE_COOKIE;
    852:     col_lens[1] = m;
    853:     col_lens[2] = n;
    854:     col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
    855:     PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
    856:     PetscFree(col_lens);

    858:     /* write out matrix, by rows */
    859:     PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);
    860:     v    = a->v;
    861:     for (i=0; i<m; i++) {
    862:       for (j=0; j<n; j++) {
    863:         vals[i + j*m] = *v++;
    864:       }
    865:     }
    866:     PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
    867:     PetscFree(vals);
    868:   } else {
    869:     PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);
    870:     col_lens[0] = MAT_FILE_COOKIE;
    871:     col_lens[1] = m;
    872:     col_lens[2] = n;
    873:     col_lens[3] = nz;

    875:     /* store lengths of each row and write (including header) to file */
    876:     for (i=0; i<m; i++) col_lens[4+i] = n;
    877:     PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);

    879:     /* Possibly should write in smaller increments, not whole matrix at once? */
    880:     /* store column indices (zero start index) */
    881:     ict = 0;
    882:     for (i=0; i<m; i++) {
    883:       for (j=0; j<n; j++) col_lens[ict++] = j;
    884:     }
    885:     PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
    886:     PetscFree(col_lens);

    888:     /* store nonzero values */
    889:     PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);
    890:     ict  = 0;
    891:     for (i=0; i<m; i++) {
    892:       v = a->v + i;
    893:       for (j=0; j<n; j++) {
    894:         anonz[ict++] = *v; v += a->lda;
    895:       }
    896:     }
    897:     PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
    898:     PetscFree(anonz);
    899:   }
    900:   return(0);
    901: }

    905: PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
    906: {
    907:   Mat               A = (Mat) Aa;
    908:   Mat_SeqDense      *a = (Mat_SeqDense*)A->data;
    909:   PetscErrorCode    ierr;
    910:   PetscInt          m = A->m,n = A->n,color,i,j;
    911:   PetscScalar       *v = a->v;
    912:   PetscViewer       viewer;
    913:   PetscDraw         popup;
    914:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0;
    915:   PetscViewerFormat format;

    919:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
    920:   PetscViewerGetFormat(viewer,&format);
    921:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

    923:   /* Loop over matrix elements drawing boxes */
    924:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
    925:     /* Blue for negative and Red for positive */
    926:     color = PETSC_DRAW_BLUE;
    927:     for(j = 0; j < n; j++) {
    928:       x_l = j;
    929:       x_r = x_l + 1.0;
    930:       for(i = 0; i < m; i++) {
    931:         y_l = m - i - 1.0;
    932:         y_r = y_l + 1.0;
    933: #if defined(PETSC_USE_COMPLEX)
    934:         if (PetscRealPart(v[j*m+i]) >  0.) {
    935:           color = PETSC_DRAW_RED;
    936:         } else if (PetscRealPart(v[j*m+i]) <  0.) {
    937:           color = PETSC_DRAW_BLUE;
    938:         } else {
    939:           continue;
    940:         }
    941: #else
    942:         if (v[j*m+i] >  0.) {
    943:           color = PETSC_DRAW_RED;
    944:         } else if (v[j*m+i] <  0.) {
    945:           color = PETSC_DRAW_BLUE;
    946:         } else {
    947:           continue;
    948:         }
    949: #endif
    950:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
    951:       }
    952:     }
    953:   } else {
    954:     /* use contour shading to indicate magnitude of values */
    955:     /* first determine max of all nonzero values */
    956:     for(i = 0; i < m*n; i++) {
    957:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
    958:     }
    959:     scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
    960:     PetscDrawGetPopup(draw,&popup);
    961:     if (popup) {PetscDrawScalePopup(popup,0.0,maxv);}
    962:     for(j = 0; j < n; j++) {
    963:       x_l = j;
    964:       x_r = x_l + 1.0;
    965:       for(i = 0; i < m; i++) {
    966:         y_l   = m - i - 1.0;
    967:         y_r   = y_l + 1.0;
    968:         color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i]));
    969:         PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
    970:       }
    971:     }
    972:   }
    973:   return(0);
    974: }

    978: PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
    979: {
    980:   PetscDraw      draw;
    981:   PetscTruth     isnull;
    982:   PetscReal      xr,yr,xl,yl,h,w;

    986:   PetscViewerDrawGetDraw(viewer,0,&draw);
    987:   PetscDrawIsNull(draw,&isnull);
    988:   if (isnull) return(0);

    990:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
    991:   xr  = A->n; yr = A->m; h = yr/10.0; w = xr/10.0;
    992:   xr += w;    yr += h;  xl = -w;     yl = -h;
    993:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
    994:   PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
    995:   PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);
    996:   return(0);
    997: }

    1001: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
    1002: {
    1004:   PetscTruth     issocket,iascii,isbinary,isdraw;

    1007:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);
    1008:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    1009:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
    1010:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);

    1012:   if (iascii) {
    1013:     MatView_SeqDense_ASCII(A,viewer);
    1014: #if defined(PETSC_USE_SOCKET_VIEWER)
    1015:   } else if (issocket) {
    1016:     Mat_SeqDense   *a = (Mat_SeqDense*)A->data;
    1017:     if (a->lda>A->m) SETERRQ(PETSC_ERR_SUP,"Case can not handle LDA");
    1018:     PetscViewerSocketPutScalar(viewer,A->m,A->n,a->v);
    1019: #endif
    1020:   } else if (isbinary) {
    1021:     MatView_SeqDense_Binary(A,viewer);
    1022:   } else if (isdraw) {
    1023:     MatView_SeqDense_Draw(A,viewer);
    1024:   } else {
    1025:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name);
    1026:   }
    1027:   return(0);
    1028: }

    1032: PetscErrorCode MatDestroy_SeqDense(Mat mat)
    1033: {
    1034:   Mat_SeqDense   *l = (Mat_SeqDense*)mat->data;

    1038: #if defined(PETSC_USE_LOG)
    1039:   PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->m,mat->n);
    1040: #endif
    1041:   if (l->pivots) {PetscFree(l->pivots);}
    1042:   if (!l->user_alloc) {PetscFree(l->v);}
    1043:   PetscFree(l);
    1044:   PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);
    1045:   return(0);
    1046: }

    1050: PetscErrorCode MatTranspose_SeqDense(Mat A,Mat *matout)
    1051: {
    1052:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    1054:   PetscInt       k,j,m,n,M;
    1055:   PetscScalar    *v,tmp;

    1058:   v = mat->v; m = A->m; M = mat->lda; n = A->n;
    1059:   if (!matout) { /* in place transpose */
    1060:     if (m != n) {
    1061:       SETERRQ(PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
    1062:     } else {
    1063:       for (j=0; j<m; j++) {
    1064:         for (k=0; k<j; k++) {
    1065:           tmp = v[j + k*M];
    1066:           v[j + k*M] = v[k + j*M];
    1067:           v[k + j*M] = tmp;
    1068:         }
    1069:       }
    1070:     }
    1071:   } else { /* out-of-place transpose */
    1072:     Mat          tmat;
    1073:     Mat_SeqDense *tmatd;
    1074:     PetscScalar  *v2;

    1076:     MatCreate(A->comm,&tmat);
    1077:     MatSetSizes(tmat,A->n,A->m,A->n,A->m);
    1078:     MatSetType(tmat,A->type_name);
    1079:     MatSeqDenseSetPreallocation(tmat,PETSC_NULL);
    1080:     tmatd = (Mat_SeqDense*)tmat->data;
    1081:     v = mat->v; v2 = tmatd->v;
    1082:     for (j=0; j<n; j++) {
    1083:       for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M];
    1084:     }
    1085:     MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
    1086:     MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
    1087:     *matout = tmat;
    1088:   }
    1089:   return(0);
    1090: }

    1094: PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg)
    1095: {
    1096:   Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
    1097:   Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
    1098:   PetscInt     i,j;
    1099:   PetscScalar  *v1 = mat1->v,*v2 = mat2->v;

    1102:   if (A1->m != A2->m) {*flg = PETSC_FALSE; return(0);}
    1103:   if (A1->n != A2->n) {*flg = PETSC_FALSE; return(0);}
    1104:   for (i=0; i<A1->m; i++) {
    1105:     v1 = mat1->v+i; v2 = mat2->v+i;
    1106:     for (j=0; j<A1->n; j++) {
    1107:       if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
    1108:       v1 += mat1->lda; v2 += mat2->lda;
    1109:     }
    1110:   }
    1111:   *flg = PETSC_TRUE;
    1112:   return(0);
    1113: }

    1117: PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
    1118: {
    1119:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    1121:   PetscInt       i,n,len;
    1122:   PetscScalar    *x,zero = 0.0;

    1125:   VecSet(v,zero);
    1126:   VecGetSize(v,&n);
    1127:   VecGetArray(v,&x);
    1128:   len = PetscMin(A->m,A->n);
    1129:   if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
    1130:   for (i=0; i<len; i++) {
    1131:     x[i] = mat->v[i*mat->lda + i];
    1132:   }
    1133:   VecRestoreArray(v,&x);
    1134:   return(0);
    1135: }

    1139: PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
    1140: {
    1141:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    1142:   PetscScalar    *l,*r,x,*v;
    1144:   PetscInt       i,j,m = A->m,n = A->n;

    1147:   if (ll) {
    1148:     VecGetSize(ll,&m);
    1149:     VecGetArray(ll,&l);
    1150:     if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
    1151:     for (i=0; i<m; i++) {
    1152:       x = l[i];
    1153:       v = mat->v + i;
    1154:       for (j=0; j<n; j++) { (*v) *= x; v+= m;}
    1155:     }
    1156:     VecRestoreArray(ll,&l);
    1157:     PetscLogFlops(n*m);
    1158:   }
    1159:   if (rr) {
    1160:     VecGetSize(rr,&n);
    1161:     VecGetArray(rr,&r);
    1162:     if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
    1163:     for (i=0; i<n; i++) {
    1164:       x = r[i];
    1165:       v = mat->v + i*m;
    1166:       for (j=0; j<m; j++) { (*v++) *= x;}
    1167:     }
    1168:     VecRestoreArray(rr,&r);
    1169:     PetscLogFlops(n*m);
    1170:   }
    1171:   return(0);
    1172: }

    1176: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
    1177: {
    1178:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
    1179:   PetscScalar  *v = mat->v;
    1180:   PetscReal    sum = 0.0;
    1181:   PetscInt     lda=mat->lda,m=A->m,i,j;

    1185:   if (type == NORM_FROBENIUS) {
    1186:     if (lda>m) {
    1187:       for (j=0; j<A->n; j++) {
    1188:         v = mat->v+j*lda;
    1189:         for (i=0; i<m; i++) {
    1190: #if defined(PETSC_USE_COMPLEX)
    1191:           sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
    1192: #else
    1193:           sum += (*v)*(*v); v++;
    1194: #endif
    1195:         }
    1196:       }
    1197:     } else {
    1198:       for (i=0; i<A->n*A->m; i++) {
    1199: #if defined(PETSC_USE_COMPLEX)
    1200:         sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
    1201: #else
    1202:         sum += (*v)*(*v); v++;
    1203: #endif
    1204:       }
    1205:     }
    1206:     *nrm = sqrt(sum);
    1207:     PetscLogFlops(2*A->n*A->m);
    1208:   } else if (type == NORM_1) {
    1209:     *nrm = 0.0;
    1210:     for (j=0; j<A->n; j++) {
    1211:       v = mat->v + j*mat->lda;
    1212:       sum = 0.0;
    1213:       for (i=0; i<A->m; i++) {
    1214:         sum += PetscAbsScalar(*v);  v++;
    1215:       }
    1216:       if (sum > *nrm) *nrm = sum;
    1217:     }
    1218:     PetscLogFlops(A->n*A->m);
    1219:   } else if (type == NORM_INFINITY) {
    1220:     *nrm = 0.0;
    1221:     for (j=0; j<A->m; j++) {
    1222:       v = mat->v + j;
    1223:       sum = 0.0;
    1224:       for (i=0; i<A->n; i++) {
    1225:         sum += PetscAbsScalar(*v); v += mat->lda;
    1226:       }
    1227:       if (sum > *nrm) *nrm = sum;
    1228:     }
    1229:     PetscLogFlops(A->n*A->m);
    1230:   } else {
    1231:     SETERRQ(PETSC_ERR_SUP,"No two norm");
    1232:   }
    1233:   return(0);
    1234: }

    1238: PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op)
    1239: {
    1240:   Mat_SeqDense   *aij = (Mat_SeqDense*)A->data;
    1242:
    1244:   switch (op) {
    1245:   case MAT_ROW_ORIENTED:
    1246:     aij->roworiented = PETSC_TRUE;
    1247:     break;
    1248:   case MAT_COLUMN_ORIENTED:
    1249:     aij->roworiented = PETSC_FALSE;
    1250:     break;
    1251:   case MAT_ROWS_SORTED:
    1252:   case MAT_ROWS_UNSORTED:
    1253:   case MAT_COLUMNS_SORTED:
    1254:   case MAT_COLUMNS_UNSORTED:
    1255:   case MAT_NO_NEW_NONZERO_LOCATIONS:
    1256:   case MAT_YES_NEW_NONZERO_LOCATIONS:
    1257:   case MAT_NEW_NONZERO_LOCATION_ERR:
    1258:   case MAT_NO_NEW_DIAGONALS:
    1259:   case MAT_YES_NEW_DIAGONALS:
    1260:   case MAT_IGNORE_OFF_PROC_ENTRIES:
    1261:   case MAT_USE_HASH_TABLE:
    1262:     PetscLogInfo((A,"MatSetOption_SeqDense:Option ignored\n"));
    1263:     break;
    1264:   case MAT_SYMMETRIC:
    1265:   case MAT_STRUCTURALLY_SYMMETRIC:
    1266:   case MAT_NOT_SYMMETRIC:
    1267:   case MAT_NOT_STRUCTURALLY_SYMMETRIC:
    1268:   case MAT_HERMITIAN:
    1269:   case MAT_NOT_HERMITIAN:
    1270:   case MAT_SYMMETRY_ETERNAL:
    1271:   case MAT_NOT_SYMMETRY_ETERNAL:
    1272:     break;
    1273:   default:
    1274:     SETERRQ(PETSC_ERR_SUP,"unknown option");
    1275:   }
    1276:   return(0);
    1277: }

    1281: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
    1282: {
    1283:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
    1285:   PetscInt       lda=l->lda,m=A->m,j;

    1288:   if (lda>m) {
    1289:     for (j=0; j<A->n; j++) {
    1290:       PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
    1291:     }
    1292:   } else {
    1293:     PetscMemzero(l->v,A->m*A->n*sizeof(PetscScalar));
    1294:   }
    1295:   return(0);
    1296: }

    1300: PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag)
    1301: {
    1302:   Mat_SeqDense   *l = (Mat_SeqDense*)A->data;
    1303:   PetscInt       n = A->n,i,j;
    1304:   PetscScalar    *slot;

    1307:   for (i=0; i<N; i++) {
    1308:     slot = l->v + rows[i];
    1309:     for (j=0; j<n; j++) { *slot = 0.0; slot += n;}
    1310:   }
    1311:   if (diag != 0.0) {
    1312:     for (i=0; i<N; i++) {
    1313:       slot = l->v + (n+1)*rows[i];
    1314:       *slot = diag;
    1315:     }
    1316:   }
    1317:   return(0);
    1318: }

    1322: PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[])
    1323: {
    1324:   Mat_SeqDense *mat = (Mat_SeqDense*)A->data;

    1327:   *array = mat->v;
    1328:   return(0);
    1329: }

    1333: PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[])
    1334: {
    1336:   *array = 0; /* user cannot accidently use the array later */
    1337:   return(0);
    1338: }

    1342: static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
    1343: {
    1344:   Mat_SeqDense   *mat = (Mat_SeqDense*)A->data;
    1346:   PetscInt       i,j,m = A->m,*irow,*icol,nrows,ncols;
    1347:   PetscScalar    *av,*bv,*v = mat->v;
    1348:   Mat            newmat;

    1351:   ISGetIndices(isrow,&irow);
    1352:   ISGetIndices(iscol,&icol);
    1353:   ISGetLocalSize(isrow,&nrows);
    1354:   ISGetLocalSize(iscol,&ncols);
    1355:
    1356:   /* Check submatrixcall */
    1357:   if (scall == MAT_REUSE_MATRIX) {
    1358:     PetscInt n_cols,n_rows;
    1359:     MatGetSize(*B,&n_rows,&n_cols);
    1360:     if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size");
    1361:     newmat = *B;
    1362:   } else {
    1363:     /* Create and fill new matrix */
    1364:     MatCreate(A->comm,&newmat);
    1365:     MatSetSizes(newmat,nrows,ncols,nrows,ncols);
    1366:     MatSetType(newmat,A->type_name);
    1367:     MatSeqDenseSetPreallocation(newmat,PETSC_NULL);
    1368:   }

    1370:   /* Now extract the data pointers and do the copy,column at a time */
    1371:   bv = ((Mat_SeqDense*)newmat->data)->v;
    1372:
    1373:   for (i=0; i<ncols; i++) {
    1374:     av = v + m*icol[i];
    1375:     for (j=0; j<nrows; j++) {
    1376:       *bv++ = av[irow[j]];
    1377:     }
    1378:   }

    1380:   /* Assemble the matrices so that the correct flags are set */
    1381:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
    1382:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);

    1384:   /* Free work space */
    1385:   ISRestoreIndices(isrow,&irow);
    1386:   ISRestoreIndices(iscol,&icol);
    1387:   *B = newmat;
    1388:   return(0);
    1389: }

    1393: PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
    1394: {
    1396:   PetscInt       i;

    1399:   if (scall == MAT_INITIAL_MATRIX) {
    1400:     PetscMalloc((n+1)*sizeof(Mat),B);
    1401:   }

    1403:   for (i=0; i<n; i++) {
    1404:     MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
    1405:   }
    1406:   return(0);
    1407: }

    1411: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
    1412: {
    1413:   Mat_SeqDense   *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data;
    1415:   PetscInt       lda1=a->lda,lda2=b->lda, m=A->m,n=A->n, j;

    1418:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
    1419:   if (A->ops->copy != B->ops->copy) {
    1420:     MatCopy_Basic(A,B,str);
    1421:     return(0);
    1422:   }
    1423:   if (m != B->m || n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
    1424:   if (lda1>m || lda2>m) {
    1425:     for (j=0; j<n; j++) {
    1426:       PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
    1427:     }
    1428:   } else {
    1429:     PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(PetscScalar));
    1430:   }
    1431:   return(0);
    1432: }

    1436: PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A)
    1437: {

    1441:    MatSeqDenseSetPreallocation(A,0);
    1442:   return(0);
    1443: }

    1445: /* -------------------------------------------------------------------*/
    1446: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
    1447:        MatGetRow_SeqDense,
    1448:        MatRestoreRow_SeqDense,
    1449:        MatMult_SeqDense,
    1450: /* 4*/ MatMultAdd_SeqDense,
    1451:        MatMultTranspose_SeqDense,
    1452:        MatMultTransposeAdd_SeqDense,
    1453:        MatSolve_SeqDense,
    1454:        MatSolveAdd_SeqDense,
    1455:        MatSolveTranspose_SeqDense,
    1456: /*10*/ MatSolveTransposeAdd_SeqDense,
    1457:        MatLUFactor_SeqDense,
    1458:        MatCholeskyFactor_SeqDense,
    1459:        MatRelax_SeqDense,
    1460:        MatTranspose_SeqDense,
    1461: /*15*/ MatGetInfo_SeqDense,
    1462:        MatEqual_SeqDense,
    1463:        MatGetDiagonal_SeqDense,
    1464:        MatDiagonalScale_SeqDense,
    1465:        MatNorm_SeqDense,
    1466: /*20*/ 0,
    1467:        0,
    1468:        0,
    1469:        MatSetOption_SeqDense,
    1470:        MatZeroEntries_SeqDense,
    1471: /*25*/ MatZeroRows_SeqDense,
    1472:        MatLUFactorSymbolic_SeqDense,
    1473:        MatLUFactorNumeric_SeqDense,
    1474:        MatCholeskyFactorSymbolic_SeqDense,
    1475:        MatCholeskyFactorNumeric_SeqDense,
    1476: /*30*/ MatSetUpPreallocation_SeqDense,
    1477:        0,
    1478:        0,
    1479:        MatGetArray_SeqDense,
    1480:        MatRestoreArray_SeqDense,
    1481: /*35*/ MatDuplicate_SeqDense,
    1482:        0,
    1483:        0,
    1484:        0,
    1485:        0,
    1486: /*40*/ MatAXPY_SeqDense,
    1487:        MatGetSubMatrices_SeqDense,
    1488:        0,
    1489:        MatGetValues_SeqDense,
    1490:        MatCopy_SeqDense,
    1491: /*45*/ 0,
    1492:        MatScale_SeqDense,
    1493:        0,
    1494:        0,
    1495:        0,
    1496: /*50*/ 0,
    1497:        0,
    1498:        0,
    1499:        0,
    1500:        0,
    1501: /*55*/ 0,
    1502:        0,
    1503:        0,
    1504:        0,
    1505:        0,
    1506: /*60*/ 0,
    1507:        MatDestroy_SeqDense,
    1508:        MatView_SeqDense,
    1509:        MatGetPetscMaps_Petsc,
    1510:        0,
    1511: /*65*/ 0,
    1512:        0,
    1513:        0,
    1514:        0,
    1515:        0,
    1516: /*70*/ 0,
    1517:        0,
    1518:        0,
    1519:        0,
    1520:        0,
    1521: /*75*/ 0,
    1522:        0,
    1523:        0,
    1524:        0,
    1525:        0,
    1526: /*80*/ 0,
    1527:        0,
    1528:        0,
    1529:        0,
    1530: /*84*/ MatLoad_SeqDense,
    1531:        0,
    1532:        0,
    1533:        0,
    1534:        0,
    1535:        0,
    1536: /*90*/ 0,
    1537:        0,
    1538:        0,
    1539:        0,
    1540:        0,
    1541: /*95*/ 0,
    1542:        0,
    1543:        0,
    1544:        0};

    1548: /*@C
    1549:    MatCreateSeqDense - Creates a sequential dense matrix that
    1550:    is stored in column major order (the usual Fortran 77 manner). Many
    1551:    of the matrix operations use the BLAS and LAPACK routines.

    1553:    Collective on MPI_Comm

    1555:    Input Parameters:
    1556: +  comm - MPI communicator, set to PETSC_COMM_SELF
    1557: .  m - number of rows
    1558: .  n - number of columns
    1559: -  data - optional location of matrix data.  Set data=PETSC_NULL for PETSc
    1560:    to control all matrix memory allocation.

    1562:    Output Parameter:
    1563: .  A - the matrix

    1565:    Notes:
    1566:    The data input variable is intended primarily for Fortran programmers
    1567:    who wish to allocate their own matrix memory space.  Most users should
    1568:    set data=PETSC_NULL.

    1570:    Level: intermediate

    1572: .keywords: dense, matrix, LAPACK, BLAS

    1574: .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues()
    1575: @*/
    1576: PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
    1577: {

    1581:   MatCreate(comm,A);
    1582:   MatSetSizes(*A,m,n,m,n);
    1583:   MatSetType(*A,MATSEQDENSE);
    1584:   MatSeqDenseSetPreallocation(*A,data);
    1585:   return(0);
    1586: }

    1590: /*@C
    1591:    MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements

    1593:    Collective on MPI_Comm

    1595:    Input Parameters:
    1596: +  A - the matrix
    1597: -  data - the array (or PETSC_NULL)

    1599:    Notes:
    1600:    The data input variable is intended primarily for Fortran programmers
    1601:    who wish to allocate their own matrix memory space.  Most users should
    1602:    set data=PETSC_NULL.

    1604:    Level: intermediate

    1606: .keywords: dense, matrix, LAPACK, BLAS

    1608: .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues()
    1609: @*/
    1610: PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
    1611: {
    1612:   PetscErrorCode ierr,(*f)(Mat,PetscScalar[]);

    1615:   PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);
    1616:   if (f) {
    1617:     (*f)(B,data);
    1618:   }
    1619:   return(0);
    1620: }

    1625: PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
    1626: {
    1627:   Mat_SeqDense   *b;

    1631:   B->preallocated = PETSC_TRUE;
    1632:   b               = (Mat_SeqDense*)B->data;
    1633:   if (!data) {
    1634:     PetscMalloc((B->m*B->n+1)*sizeof(PetscScalar),&b->v);
    1635:     PetscMemzero(b->v,B->m*B->n*sizeof(PetscScalar));
    1636:     b->user_alloc = PETSC_FALSE;
    1637:     PetscLogObjectMemory(B,B->n*B->m*sizeof(PetscScalar));
    1638:   } else { /* user-allocated storage */
    1639:     b->v          = data;
    1640:     b->user_alloc = PETSC_TRUE;
    1641:   }
    1642:   return(0);
    1643: }

    1648: /*@C
    1649:   MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array

    1651:   Input parameter:
    1652: + A - the matrix
    1653: - lda - the leading dimension

    1655:   Notes:
    1656:   This routine is to be used in conjunction with MatSeqDenseSetPreallocation;
    1657:   it asserts that the preallocation has a leading dimension (the LDA parameter
    1658:   of Blas and Lapack fame) larger than M, the first dimension of the matrix.

    1660:   Level: intermediate

    1662: .keywords: dense, matrix, LAPACK, BLAS

    1664: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation()
    1665: @*/
    1666: PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda)
    1667: {
    1668:   Mat_SeqDense *b = (Mat_SeqDense*)B->data;
    1670:   if (lda < B->m) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->m);
    1671:   b->lda = lda;
    1672:   return(0);
    1673: }

    1675: /*MC
    1676:    MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.

    1678:    Options Database Keys:
    1679: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()

    1681:   Level: beginner

    1683: .seealso: MatCreateSeqDense
    1684: M*/

    1689: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B)
    1690: {
    1691:   Mat_SeqDense   *b;
    1693:   PetscMPIInt    size;

    1696:   MPI_Comm_size(B->comm,&size);
    1697:   if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1");

    1699:   B->m = B->M = PetscMax(B->m,B->M);
    1700:   B->n = B->N = PetscMax(B->n,B->N);

    1702:   PetscNew(Mat_SeqDense,&b);
    1703:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
    1704:   B->factor       = 0;
    1705:   B->mapping      = 0;
    1706:   PetscLogObjectMemory(B,sizeof(struct _p_Mat));
    1707:   B->data         = (void*)b;

    1709:   PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);
    1710:   PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);

    1712:   b->pivots       = 0;
    1713:   b->roworiented  = PETSC_TRUE;
    1714:   b->v            = 0;
    1715:   b->lda          = B->m;

    1717:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C",
    1718:                                     "MatSeqDenseSetPreallocation_SeqDense",
    1719:                                      MatSeqDenseSetPreallocation_SeqDense);
    1720:   return(0);
    1721: }