Actual source code: aijsbaij.c

      2: #include <../src/mat/impls/aij/seq/aij.h>
      3: #include <../src/mat/impls/baij/seq/baij.h>
      4: #include <../src/mat/impls/sbaij/seq/sbaij.h>

      6: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
      7: {
      8:   Mat            B;
      9:   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data;
     10:   Mat_SeqAIJ     *b;
     12:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp;
     13:   PetscInt       bs =A->rmap->bs,bs2=bs*bs,mbs=A->rmap->N/bs,diagcnt=0;
     14:   MatScalar      *av,*bv;
     15: #if defined(PETSC_USE_COMPLEX)
     16:   const int      aconj = A->hermitian ? 1 : 0;
     17: #else
     18:   const int      aconj = 0;
     19: #endif

     22:   /* compute rowlengths of newmat */
     23:   PetscMalloc2(m,&rowlengths,m+1,&rowstart);

     25:   for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
     26:   k  = 0;
     27:   for (i=0; i<mbs; i++) {
     28:     nz = ai[i+1] - ai[i];
     29:     if (nz) {
     30:       rowlengths[k] += nz;   /* no. of upper triangular blocks */
     31:       if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */
     32:       for (j=0; j<nz; j++) { /* no. of lower triangular blocks */
     33:         rowlengths[(*aj)*bs]++; aj++;
     34:       }
     35:     }
     36:     rowlengths[k] *= bs;
     37:     for (j=1; j<bs; j++) {
     38:       rowlengths[k+j] = rowlengths[k];
     39:     }
     40:     k += bs;
     41:   }

     43:   if (reuse != MAT_REUSE_MATRIX) {
     44:     MatCreate(PetscObjectComm((PetscObject)A),&B);
     45:     MatSetSizes(B,m,n,m,n);
     46:     MatSetType(B,MATSEQAIJ);
     47:     MatSeqAIJSetPreallocation(B,0,rowlengths);
     48:     MatSetBlockSize(B,A->rmap->bs);
     49:   } else B = *newmat;

     51:   b  = (Mat_SeqAIJ*)(B->data);
     52:   bi = b->i;
     53:   bj = b->j;
     54:   bv = b->a;

     56:   /* set b->i */
     57:   bi[0] = 0; rowstart[0] = 0;
     58:   for (i=0; i<mbs; i++) {
     59:     for (j=0; j<bs; j++) {
     60:       b->ilen[i*bs+j]    = rowlengths[i*bs];
     61:       rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs];
     62:     }
     63:     bi[i+1] = bi[i] + rowlengths[i*bs]/bs;
     64:   }
     65:   if (bi[mbs] != 2*a->nz - diagcnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz-diagcnt: %D\n",bi[mbs],2*a->nz - diagcnt);

     67:   /* set b->j and b->a */
     68:   aj = a->j; av = a->a;
     69:   for (i=0; i<mbs; i++) {
     70:     nz = ai[i+1] - ai[i];
     71:     /* diagonal block */
     72:     if (nz && *aj == i) {
     73:       nz--;
     74:       for (j=0; j<bs; j++) {   /* row i*bs+j */
     75:         itmp = i*bs+j;
     76:         for (k=0; k<bs; k++) { /* col i*bs+k */
     77:           *(bj + rowstart[itmp]) = (*aj)*bs+k;
     78:           *(bv + rowstart[itmp]) = *(av+k*bs+j);
     79:           rowstart[itmp]++;
     80:         }
     81:       }
     82:       aj++; av += bs2;
     83:     }

     85:     while (nz--) {
     86:       /* lower triangular blocks */
     87:       for (j=0; j<bs; j++) {   /* row (*aj)*bs+j */
     88:         itmp = (*aj)*bs+j;
     89:         for (k=0; k<bs; k++) { /* col i*bs+k */
     90:           *(bj + rowstart[itmp]) = i*bs+k;
     91:           *(bv + rowstart[itmp]) = aconj ? PetscConj(*(av+j*bs+k)) : *(av+j*bs+k);
     92:           rowstart[itmp]++;
     93:         }
     94:       }
     95:       /* upper triangular blocks */
     96:       for (j=0; j<bs; j++) {   /* row i*bs+j */
     97:         itmp = i*bs+j;
     98:         for (k=0; k<bs; k++) { /* col (*aj)*bs+k */
     99:           *(bj + rowstart[itmp]) = (*aj)*bs+k;
    100:           *(bv + rowstart[itmp]) = *(av+k*bs+j);
    101:           rowstart[itmp]++;
    102:         }
    103:       }
    104:       aj++; av += bs2;
    105:     }
    106:   }
    107:   PetscFree2(rowlengths,rowstart);
    108:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    109:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

    111:   if (reuse == MAT_INPLACE_MATRIX) {
    112:     MatHeaderReplace(A,&B);
    113:   } else {
    114:     *newmat = B;
    115:   }
    116:   return(0);
    117: }

    119: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
    120: {
    121:   Mat            B;
    122:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
    123:   Mat_SeqSBAIJ   *b;
    125:   PetscInt       *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths,bs=PetscAbs(A->rmap->bs);
    126:   MatScalar      *av,*bv;
    127:   PetscBool      miss = PETSC_FALSE;

    130: #if !defined(PETSC_USE_COMPLEX)
    131:   if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
    132: #else
    133:   if (!A->symmetric && !A->hermitian) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be either symmetric or hermitian. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE) and/or MatSetOption(mat,MAT_HERMITIAN,PETSC_TRUE)");
    134: #endif
    135:   if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");

    137:   PetscMalloc1(m/bs,&rowlengths);
    138:   for (i=0; i<m/bs; i++) {
    139:     if (a->diag[i*bs] == ai[i*bs+1]) { /* missing diagonal */
    140:       rowlengths[i] = (ai[i*bs+1] - ai[i*bs])/bs; /* allocate some extra space */
    141:       miss = PETSC_TRUE;
    142:     } else {
    143:       rowlengths[i] = (ai[i*bs+1] - a->diag[i*bs])/bs;
    144:     }
    145:   }
    146:   if (reuse != MAT_REUSE_MATRIX) {
    147:     MatCreate(PetscObjectComm((PetscObject)A),&B);
    148:     MatSetSizes(B,m,n,m,n);
    149:     MatSetType(B,MATSEQSBAIJ);
    150:     MatSeqSBAIJSetPreallocation(B,bs,0,rowlengths);
    151:   } else B = *newmat;

    153:   if (bs == 1 && !miss) {
    154:     b  = (Mat_SeqSBAIJ*)(B->data);
    155:     bi = b->i;
    156:     bj = b->j;
    157:     bv = b->a;

    159:     bi[0] = 0;
    160:     for (i=0; i<m; i++) {
    161:       aj = a->j + a->diag[i];
    162:       av = a->a + a->diag[i];
    163:       for (j=0; j<rowlengths[i]; j++) {
    164:         *bj = *aj; bj++; aj++;
    165:         *bv = *av; bv++; av++;
    166:       }
    167:       bi[i+1]    = bi[i] + rowlengths[i];
    168:       b->ilen[i] = rowlengths[i];
    169:     }
    170:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    171:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    172:   } else {
    173:     MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
    174:     /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
    175:     /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before                                          */
    176:     /*                      MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below                                        */
    177:     MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);
    178:   }
    179:   PetscFree(rowlengths);
    180:   if (reuse == MAT_INPLACE_MATRIX) {
    181:     MatHeaderReplace(A,&B);
    182:   } else *newmat = B;
    183:   return(0);
    184: }

    186: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
    187: {
    188:   Mat            B;
    189:   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data;
    190:   Mat_SeqBAIJ    *b;
    192:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
    193:   PetscInt       bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
    194:   MatScalar      *av,*bv;
    195: #if defined(PETSC_USE_COMPLEX)
    196:   const int      aconj = A->hermitian ? 1 : 0;
    197: #else
    198:   const int      aconj = 0;
    199: #endif

    202:   /* compute browlengths of newmat */
    203:   PetscMalloc2(mbs,&browlengths,mbs,&browstart);
    204:   for (i=0; i<mbs; i++) browlengths[i] = 0;
    205:   for (i=0; i<mbs; i++) {
    206:     nz = ai[i+1] - ai[i];
    207:     aj++; /* skip diagonal */
    208:     for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
    209:       browlengths[*aj]++; aj++;
    210:     }
    211:     browlengths[i] += nz;   /* no. of upper triangular blocks */
    212:   }

    214:   if (reuse != MAT_REUSE_MATRIX) {
    215:     MatCreate(PetscObjectComm((PetscObject)A),&B);
    216:     MatSetSizes(B,m,n,m,n);
    217:     MatSetType(B,MATSEQBAIJ);
    218:     MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
    219:   } else B = *newmat;

    221:   b  = (Mat_SeqBAIJ*)(B->data);
    222:   bi = b->i;
    223:   bj = b->j;
    224:   bv = b->a;

    226:   /* set b->i */
    227:   bi[0] = 0;
    228:   for (i=0; i<mbs; i++) {
    229:     b->ilen[i]   = browlengths[i];
    230:     bi[i+1]      = bi[i] + browlengths[i];
    231:     browstart[i] = bi[i];
    232:   }
    233:   if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz - mbs: %D\n",bi[mbs],2*a->nz - mbs);

    235:   /* set b->j and b->a */
    236:   aj = a->j; av = a->a;
    237:   for (i=0; i<mbs; i++) {
    238:     /* diagonal block */
    239:     *(bj + browstart[i]) = *aj; aj++;

    241:     itmp = bs2*browstart[i];
    242:     for (k=0; k<bs2; k++) {
    243:       *(bv + itmp + k) = *av; av++;
    244:     }
    245:     browstart[i]++;

    247:     nz = ai[i+1] - ai[i] -1;
    248:     while (nz--) {
    249:       /* lower triangular blocks - transpose blocks of A */
    250:       *(bj + browstart[*aj]) = i; /* block col index */

    252:       itmp = bs2*browstart[*aj];  /* row index */
    253:       for (col=0; col<bs; col++) {
    254:         k = col;
    255:         for (row=0; row<bs; row++) {
    256:           bv[itmp + col*bs+row] = aconj ? PetscConj(av[k]) : av[k];
    257:           k+=bs;
    258:         }
    259:       }
    260:       browstart[*aj]++;

    262:       /* upper triangular blocks */
    263:       *(bj + browstart[i]) = *aj; aj++;

    265:       itmp = bs2*browstart[i];
    266:       for (k=0; k<bs2; k++) {
    267:         bv[itmp + k] = av[k];
    268:       }
    269:       av += bs2;
    270:       browstart[i]++;
    271:     }
    272:   }
    273:   PetscFree2(browlengths,browstart);
    274:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    275:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

    277:   if (reuse == MAT_INPLACE_MATRIX) {
    278:     MatHeaderReplace(A,&B);
    279:   } else *newmat = B;
    280:   return(0);
    281: }

    283: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
    284: {
    285:   Mat            B;
    286:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
    287:   Mat_SeqSBAIJ   *b;
    289:   PetscInt       *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
    290:   PetscInt       bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
    291:   MatScalar      *av,*bv;
    292:   PetscBool      flg;

    295:   if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
    296:   if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
    297:   MatMissingDiagonal_SeqBAIJ(A,&flg,&dd); /* check for missing diagonals, then mark diag */
    298:   if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd);

    300:   PetscMalloc1(mbs,&browlengths);
    301:   for (i=0; i<mbs; i++) {
    302:     browlengths[i] = ai[i+1] - a->diag[i];
    303:   }

    305:   if (reuse != MAT_REUSE_MATRIX) {
    306:     MatCreate(PetscObjectComm((PetscObject)A),&B);
    307:     MatSetSizes(B,m,n,m,n);
    308:     MatSetType(B,MATSEQSBAIJ);
    309:     MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);
    310:   } else B = *newmat;

    312:   b  = (Mat_SeqSBAIJ*)(B->data);
    313:   bi = b->i;
    314:   bj = b->j;
    315:   bv = b->a;

    317:   bi[0] = 0;
    318:   for (i=0; i<mbs; i++) {
    319:     aj = a->j + a->diag[i];
    320:     av = a->a + (a->diag[i])*bs2;
    321:     for (j=0; j<browlengths[i]; j++) {
    322:       *bj = *aj; bj++; aj++;
    323:       for (k=0; k<bs2; k++) {
    324:         *bv = *av; bv++; av++;
    325:       }
    326:     }
    327:     bi[i+1]    = bi[i] + browlengths[i];
    328:     b->ilen[i] = browlengths[i];
    329:   }
    330:   PetscFree(browlengths);
    331:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    332:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

    334:   if (reuse == MAT_INPLACE_MATRIX) {
    335:     MatHeaderReplace(A,&B);
    336:   } else *newmat = B;
    337:   return(0);
    338: }