Actual source code: aijbaij.c

      1: #define PETSCMAT_DLL

     3:  #include src/mat/impls/baij/seq/baij.h

      8: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
      9: {
     10:   Mat            B;
     11:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
     13:   PetscInt       bs = A->bs,*ai = a->i,*aj = a->j,n = A->M/bs,i,j,k;
     14:   PetscInt       *rowlengths,*rows,*cols,maxlen = 0,ncols;
     15:   PetscScalar    *aa = a->a;

     18:   PetscMalloc(n*bs*sizeof(PetscInt),&rowlengths);
     19:   for (i=0; i<n; i++) {
     20:     maxlen = PetscMax(maxlen,(ai[i+1] - ai[i]));
     21:     for (j=0; j<bs; j++) {
     22:       rowlengths[i*bs+j] = bs*(ai[i+1] - ai[i]);
     23:     }
     24:   }
     25:   MatCreate(A->comm,&B);
     26:   MatSetSizes(B,A->m,A->n,A->m,A->n);
     27:   MatSetType(B,newtype);
     28:   MatSeqAIJSetPreallocation(B,0,rowlengths);
     29:   MatSetOption(B,MAT_COLUMN_ORIENTED);
     30:   MatSetOption(B,MAT_ROWS_SORTED);
     31:   MatSetOption(B,MAT_COLUMNS_SORTED);
     32:   PetscFree(rowlengths);

     34:   PetscMalloc(bs*sizeof(PetscInt),&rows);
     35:   PetscMalloc(bs*maxlen*sizeof(PetscInt),&cols);
     36:   for (i=0; i<n; i++) {
     37:     for (j=0; j<bs; j++) {
     38:       rows[j] = i*bs+j;
     39:     }
     40:     ncols = ai[i+1] - ai[i];
     41:     for (k=0; k<ncols; k++) {
     42:       for (j=0; j<bs; j++) {
     43:         cols[k*bs+j] = bs*(*aj) + j;
     44:       }
     45:       aj++;
     46:     }
     47:     MatSetValues(B,bs,rows,bs*ncols,cols,aa,INSERT_VALUES);
     48:     aa   += ncols*bs*bs;
     49:   }
     50:   PetscFree(cols);
     51:   PetscFree(rows);
     52:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
     53:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
     54:
     55:   B->bs = A->bs;

     57:   if (reuse == MAT_REUSE_MATRIX) {
     58:     MatHeaderReplace(A,B);
     59:   } else {
     60:     *newmat = B;
     61:   }
     62:   return(0);
     63: }

     66:  #include src/mat/impls/aij/seq/aij.h

     71: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
     72: {
     73:   Mat            B;
     74:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
     75:   Mat_SeqBAIJ    *b;
     77:   PetscInt       *ai=a->i,m=A->M,n=A->N,i,*rowlengths;

     80:   if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");

     82:   PetscMalloc(m*sizeof(PetscInt),&rowlengths);
     83:   for (i=0; i<m; i++) {
     84:     rowlengths[i] = ai[i+1] - ai[i];
     85:   }
     86:   MatCreate(A->comm,&B);
     87:   MatSetSizes(B,m,n,m,n);
     88:   MatSetType(B,newtype);
     89:   MatSeqBAIJSetPreallocation_SeqBAIJ(B,1,0,rowlengths);
     90:   PetscFree(rowlengths);

     92:   MatSetOption(B,MAT_ROW_ORIENTED);
     93:   MatSetOption(B,MAT_ROWS_SORTED);
     94:   MatSetOption(B,MAT_COLUMNS_SORTED);
     95:
     96:   b  = (Mat_SeqBAIJ*)(B->data);

     98:   PetscMemcpy(b->i,a->i,(m+1)*sizeof(PetscInt));
     99:   PetscMemcpy(b->ilen,a->ilen,m*sizeof(PetscInt));
    100:   PetscMemcpy(b->j,a->j,a->nz*sizeof(PetscInt));
    101:   PetscMemcpy(b->a,a->a,a->nz*sizeof(MatScalar));
    102:
    103:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    104:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

    106:   if (reuse == MAT_REUSE_MATRIX) {
    107:     MatHeaderReplace(A,B);
    108:   } else {
    109:     *newmat = B;
    110:   }
    111:   return(0);
    112: }