Actual source code: aijsbaij.c

      1: #define PETSCMAT_DLL

     3:  #include src/mat/impls/aij/seq/aij.h
     4:  #include src/mat/impls/baij/seq/baij.h
     5:  #include src/mat/impls/sbaij/seq/sbaij.h

     10: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
     11: {
     12:   Mat            B;
     13:   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data;
     14:   Mat_SeqAIJ     *b;
     16:   PetscInt       *ai=a->i,*aj=a->j,m=A->m,n=A->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp;
     17:   PetscInt       bs=A->bs,bs2=bs*bs,mbs=A->m/bs;
     18:   PetscScalar    *av,*bv;

     22:   /* compute rowlengths of newmat */
     23:   PetscMalloc((2*m+1)*sizeof(PetscInt),&rowlengths);
     24:   rowstart = rowlengths + m;
     25:
     26:   for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
     27:   aj = a->j;
     28:   k = 0;
     29:   for (i=0; i<mbs; i++) {
     30:     nz = ai[i+1] - ai[i];
     31:     aj++; /* skip diagonal */
     32:     for (j=1; j<nz; j++) { /* no. of lower triangular blocks */
     33:       rowlengths[(*aj)*bs]++; aj++;
     34:     }
     35:     rowlengths[k] += nz;   /* no. of upper triangular blocks */
     36:     rowlengths[k] *= bs;
     37:     for (j=1; j<bs; j++) {
     38:       rowlengths[k+j] = rowlengths[k];
     39:     }
     40:     k += bs;
     41:     /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */
     42:   }
     43:
     44:   MatCreate(A->comm,&B);
     45:   MatSetSizes(B,m,n,m,n);
     46:   MatSetType(B,newtype);
     47:   MatSeqAIJSetPreallocation(B,0,rowlengths);
     48:   MatSetOption(B,MAT_COLUMN_ORIENTED);
     49:   MatSetOption(B,MAT_ROWS_SORTED);
     50:   MatSetOption(B,MAT_COLUMNS_SORTED);
     51:   B->bs = A->bs;

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

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

     69:   /* set b->j and b->a */
     70:   aj = a->j; av = a->a;
     71:   for (i=0; i<mbs; i++) {
     72:     /* diagonal block */
     73:     for (j=0; j<bs; j++){   /* row i*bs+j */
     74:       itmp = i*bs+j;
     75:       for (k=0; k<bs; k++){ /* col i*bs+k */
     76:         *(bj + rowstart[itmp]) = (*aj)*bs+k;
     77:         *(bv + rowstart[itmp]) = *(av+k*bs+j);
     78:         rowstart[itmp]++;
     79:       }
     80:     }
     81:     aj++; av += bs2;
     82:
     83:     nz = ai[i+1] - ai[i] -1;
     84:     while (nz--){
     85:       /* lower triangular blocks */
     86:       for (j=0; j<bs; j++){   /* row (*aj)*bs+j */
     87:         itmp = (*aj)*bs+j;
     88:         for (k=0; k<bs; k++){ /* col i*bs+k */
     89:           *(bj + rowstart[itmp]) = i*bs+k;
     90:           *(bv + rowstart[itmp]) = *(av+k*bs+j);
     91:           rowstart[itmp]++;
     92:         }
     93:       }
     94:       /* upper triangular blocks */
     95:       for (j=0; j<bs; j++){   /* row i*bs+j */
     96:         itmp = i*bs+j;
     97:         for (k=0; k<bs; k++){ /* col (*aj)*bs+k */
     98:           *(bj + rowstart[itmp]) = (*aj)*bs+k;
     99:           *(bv + rowstart[itmp]) = *(av+k*bs+j);
    100:           rowstart[itmp]++;
    101:         }
    102:       }
    103:       aj++; av += bs2;
    104:     }
    105:   }
    106:   PetscFree(rowlengths);
    107:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    108:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

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

    122: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) {
    123:   Mat            B;
    124:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
    125:   Mat_SeqSBAIJ   *b;
    127:   PetscInt       *ai=a->i,*aj,m=A->M,n=A->N,i,j,*bi,*bj,*rowlengths;
    128:   PetscScalar    *av,*bv;

    131:   if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
    132:   MatMarkDiagonal_SeqAIJ(A); /* check for missing diagonals, then mark diag */

    134:   PetscMalloc(m*sizeof(PetscInt),&rowlengths);
    135:   for (i=0; i<m; i++) {
    136:     rowlengths[i] = ai[i+1] - a->diag[i];
    137:   }
    138:   MatCreate(A->comm,&B);
    139:   MatSetSizes(B,m,n,m,n);
    140:   MatSetType(B,newtype);
    141:   MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,1,0,rowlengths);

    143:   MatSetOption(B,MAT_ROW_ORIENTED);
    144:   MatSetOption(B,MAT_ROWS_SORTED);
    145:   MatSetOption(B,MAT_COLUMNS_SORTED);
    146:
    147:   b  = (Mat_SeqSBAIJ*)(B->data);
    148:   bi = b->i;
    149:   bj = b->j;
    150:   bv = b->a;

    152:   bi[0] = 0;
    153:   for (i=0; i<m; i++) {
    154:     aj = a->j + a->diag[i];
    155:     av = a->a + a->diag[i];
    156:     for (j=0; j<rowlengths[i]; j++){
    157:       *bj = *aj; bj++; aj++;
    158:       *bv = *av; bv++; av++;
    159:     }
    160:     bi[i+1]    = bi[i] + rowlengths[i];
    161:     b->ilen[i] = rowlengths[i];
    162:   }
    163:
    164:   PetscFree(rowlengths);
    165:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    166:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

    168:   if (reuse == MAT_REUSE_MATRIX) {
    169:     MatHeaderReplace(A,B);
    170:   } else {
    171:     *newmat = B;
    172:   }
    173:   return(0);
    174: }

    180: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
    181: {
    182:   Mat            B;
    183:   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data;
    184:   Mat_SeqBAIJ    *b;
    186:   PetscInt       *ai=a->i,*aj=a->j,m=A->m,n=A->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
    187:   PetscInt       bs=A->bs,bs2=bs*bs,mbs=m/bs;
    188:   PetscScalar    *av,*bv;

    191:   /* compute browlengths of newmat */
    192:   PetscMalloc(2*mbs*sizeof(PetscInt),&browlengths);
    193:   browstart = browlengths + mbs;
    194:   for (i=0; i<mbs; i++) browlengths[i] = 0;
    195:   aj = a->j;
    196:   for (i=0; i<mbs; i++) {
    197:     nz = ai[i+1] - ai[i];
    198:     aj++; /* skip diagonal */
    199:     for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
    200:       browlengths[*aj]++; aj++;
    201:     }
    202:     browlengths[i] += nz;   /* no. of upper triangular blocks */
    203:   }
    204:
    205:   MatCreate(A->comm,&B);
    206:   MatSetSizes(B,m,n,m,n);
    207:   MatSetType(B,newtype);
    208:   MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
    209:   MatSetOption(B,MAT_ROW_ORIENTED);
    210:   MatSetOption(B,MAT_ROWS_SORTED);
    211:   MatSetOption(B,MAT_COLUMNS_SORTED);
    212:
    213:   b  = (Mat_SeqBAIJ*)(B->data);
    214:   bi = b->i;
    215:   bj = b->j;
    216:   bv = b->a;

    218:   /* set b->i */
    219:   bi[0] = 0;
    220:   for (i=0; i<mbs; i++){
    221:     b->ilen[i]    = browlengths[i];
    222:     bi[i+1]       = bi[i] + browlengths[i];
    223:     browstart[i]  = bi[i];
    224:   }
    225:   if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz - mbs: %D\n",bi[mbs],2*a->nz - mbs);
    226:
    227:   /* set b->j and b->a */
    228:   aj = a->j; av = a->a;
    229:   for (i=0; i<mbs; i++) {
    230:     /* diagonal block */
    231:     *(bj + browstart[i]) = *aj; aj++;
    232:     itmp = bs2*browstart[i];
    233:     for (k=0; k<bs2; k++){
    234:       *(bv + itmp + k) = *av; av++;
    235:     }
    236:     browstart[i]++;
    237:
    238:     nz = ai[i+1] - ai[i] -1;
    239:     while (nz--){
    240:       /* lower triangular blocks */
    241:       *(bj + browstart[*aj]) = i;
    242:       itmp = bs2*browstart[*aj];
    243:       for (k=0; k<bs2; k++){
    244:         *(bv + itmp + k) = *(av + k);
    245:       }
    246:       browstart[*aj]++;

    248:       /* upper triangular blocks */
    249:       *(bj + browstart[i]) = *aj; aj++;
    250:       itmp = bs2*browstart[i];
    251:       for (k=0; k<bs2; k++){
    252:         *(bv + itmp + k) = *av; av++;
    253:       }
    254:       browstart[i]++;
    255:     }
    256:   }
    257:   PetscFree(browlengths);
    258:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    259:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

    261:   if (reuse == MAT_REUSE_MATRIX) {
    262:     MatHeaderReplace(A,B);
    263:   } else {
    264:     *newmat = B;
    265:   }
    266:   return(0);
    267: }

    273: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
    274: {
    275:   Mat            B;
    276:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
    277:   Mat_SeqSBAIJ   *b;
    279:   PetscInt       *ai=a->i,*aj,m=A->m,n=A->n,i,j,k,*bi,*bj,*browlengths;
    280:   PetscInt       bs=A->bs,bs2=bs*bs,mbs=m/bs;
    281:   PetscScalar    *av,*bv;

    284:   if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
    285:   MatMissingDiagonal_SeqBAIJ(A); /* check for missing diagonals, then mark diag */
    286:
    287:   PetscMalloc(mbs*sizeof(PetscInt),&browlengths);
    288:   for (i=0; i<mbs; i++) {
    289:     browlengths[i] = ai[i+1] - a->diag[i];
    290:   }

    292:   MatCreate(A->comm,&B);
    293:   MatSetSizes(B,m,n,m,n);
    294:   MatSetType(B,newtype);
    295:   MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,browlengths);
    296:   MatSetOption(B,MAT_ROW_ORIENTED);
    297:   MatSetOption(B,MAT_ROWS_SORTED);
    298:   MatSetOption(B,MAT_COLUMNS_SORTED);
    299:
    300:   b  = (Mat_SeqSBAIJ*)(B->data);
    301:   bi = b->i;
    302:   bj = b->j;
    303:   bv = b->a;

    305:   bi[0] = 0;
    306:   for (i=0; i<mbs; i++) {
    307:     aj = a->j + a->diag[i];
    308:     av = a->a + (a->diag[i])*bs2;
    309:     for (j=0; j<browlengths[i]; j++){
    310:       *bj = *aj; bj++; aj++;
    311:       for (k=0; k<bs2; k++){
    312:         *bv = *av; bv++; av++;
    313:       }
    314:     }
    315:     bi[i+1]    = bi[i] + browlengths[i];
    316:     b->ilen[i] = browlengths[i];
    317:   }
    318:   PetscFree(browlengths);
    319:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    320:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

    322:   if (reuse == MAT_REUSE_MATRIX) {
    323:     MatHeaderReplace(A,B);
    324:   } else {
    325:     *newmat = B;
    326:   }

    328:   return(0);
    329: }