Actual source code: aijbaij.c


  2: #include <../src/mat/impls/baij/seq/baij.h>

  4: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
  5: {
  6:   Mat          B;
  7:   Mat_SeqAIJ  *b;
  8:   PetscBool    roworiented;
  9:   Mat_SeqBAIJ *a  = (Mat_SeqBAIJ *)A->data;
 10:   PetscInt     bs = A->rmap->bs, *ai = a->i, *aj = a->j, n = A->rmap->N / bs, i, j, k;
 11:   PetscInt    *rowlengths, *rows, *cols, maxlen            = 0, ncols;
 12:   MatScalar   *aa = a->a;

 14:   if (reuse == MAT_REUSE_MATRIX) {
 15:     B = *newmat;
 16:     for (i = 0; i < n; i++) maxlen = PetscMax(maxlen, (ai[i + 1] - ai[i]));
 17:   } else {
 18:     PetscMalloc1(n * bs, &rowlengths);
 19:     for (i = 0; i < n; i++) {
 20:       maxlen = PetscMax(maxlen, (ai[i + 1] - ai[i]));
 21:       for (j = 0; j < bs; j++) rowlengths[i * bs + j] = bs * (ai[i + 1] - ai[i]);
 22:     }
 23:     MatCreate(PetscObjectComm((PetscObject)A), &B);
 24:     MatSetType(B, MATSEQAIJ);
 25:     MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
 26:     MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs);
 27:     MatSeqAIJSetPreallocation(B, 0, rowlengths);
 28:     PetscFree(rowlengths);
 29:   }
 30:   b           = (Mat_SeqAIJ *)B->data;
 31:   roworiented = b->roworiented;

 33:   MatSetOption(B, MAT_ROW_ORIENTED, PETSC_FALSE);
 34:   PetscMalloc1(bs, &rows);
 35:   PetscMalloc1(bs * maxlen, &cols);
 36:   for (i = 0; i < n; i++) {
 37:     for (j = 0; j < bs; j++) rows[j] = i * bs + j;
 38:     ncols = ai[i + 1] - ai[i];
 39:     for (k = 0; k < ncols; k++) {
 40:       for (j = 0; j < bs; j++) cols[k * bs + j] = bs * (*aj) + j;
 41:       aj++;
 42:     }
 43:     MatSetValues(B, bs, rows, bs * ncols, cols, aa, INSERT_VALUES);
 44:     aa += ncols * bs * bs;
 45:   }
 46:   PetscFree(cols);
 47:   PetscFree(rows);
 48:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 49:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 50:   MatSetOption(B, MAT_ROW_ORIENTED, roworiented);

 52:   if (reuse == MAT_INPLACE_MATRIX) {
 53:     MatHeaderReplace(A, &B);
 54:   } else *newmat = B;
 55:   return 0;
 56: }

 58: #include <../src/mat/impls/aij/seq/aij.h>

 60: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
 61: {
 62:   Mat          B;
 63:   Mat_SeqAIJ  *a = (Mat_SeqAIJ *)A->data;
 64:   Mat_SeqBAIJ *b;
 65:   PetscInt    *ai = a->i, m = A->rmap->N, n = A->cmap->N, i, *rowlengths, bs = PetscAbs(A->rmap->bs);

 67:   if (reuse != MAT_REUSE_MATRIX) {
 68:     PetscMalloc1(m / bs, &rowlengths);
 69:     for (i = 0; i < m / bs; i++) rowlengths[i] = (ai[i * bs + 1] - ai[i * bs]) / bs;
 70:     MatCreate(PetscObjectComm((PetscObject)A), &B);
 71:     MatSetSizes(B, m, n, m, n);
 72:     MatSetType(B, MATSEQBAIJ);
 73:     MatSeqBAIJSetPreallocation(B, bs, 0, rowlengths);
 74:     PetscFree(rowlengths);
 75:   } else B = *newmat;

 77:   if (bs == 1) {
 78:     b = (Mat_SeqBAIJ *)(B->data);

 80:     PetscArraycpy(b->i, a->i, m + 1);
 81:     PetscArraycpy(b->ilen, a->ilen, m);
 82:     PetscArraycpy(b->j, a->j, a->nz);
 83:     PetscArraycpy(b->a, a->a, a->nz);

 85:     MatSetOption(B, MAT_ROW_ORIENTED, PETSC_TRUE);
 86:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 87:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 88:   } else {
 89:     /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
 90:     /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before                                          */
 91:     /*                      MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below                                        */
 92:     MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B);
 93:   }

 95:   if (reuse == MAT_INPLACE_MATRIX) {
 96:     MatHeaderReplace(A, &B);
 97:   } else *newmat = B;
 98:   return 0;
 99: }