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++) {
17: maxlen = PetscMax(maxlen,(ai[i+1] - ai[i]));
18: }
19: } else {
20: PetscMalloc1(n*bs,&rowlengths);
21: for (i=0; i<n; i++) {
22: maxlen = PetscMax(maxlen,(ai[i+1] - ai[i]));
23: for (j=0; j<bs; j++) {
24: rowlengths[i*bs+j] = bs*(ai[i+1] - ai[i]);
25: }
26: }
27: MatCreate(PetscObjectComm((PetscObject)A),&B);
28: MatSetType(B,MATSEQAIJ);
29: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
30: MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
31: MatSeqAIJSetPreallocation(B,0,rowlengths);
32: PetscFree(rowlengths);
33: }
34: b = (Mat_SeqAIJ*)B->data;
35: roworiented = b->roworiented;
37: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);
38: PetscMalloc1(bs,&rows);
39: PetscMalloc1(bs*maxlen,&cols);
40: for (i=0; i<n; i++) {
41: for (j=0; j<bs; j++) {
42: rows[j] = i*bs+j;
43: }
44: ncols = ai[i+1] - ai[i];
45: for (k=0; k<ncols; k++) {
46: for (j=0; j<bs; j++) {
47: cols[k*bs+j] = bs*(*aj) + j;
48: }
49: aj++;
50: }
51: MatSetValues(B,bs,rows,bs*ncols,cols,aa,INSERT_VALUES);
52: aa += ncols*bs*bs;
53: }
54: PetscFree(cols);
55: PetscFree(rows);
56: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
57: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
58: MatSetOption(B,MAT_ROW_ORIENTED,roworiented);
60: if (reuse == MAT_INPLACE_MATRIX) {
61: MatHeaderReplace(A,&B);
62: } else *newmat = B;
63: return 0;
64: }
66: #include <../src/mat/impls/aij/seq/aij.h>
68: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
69: {
70: Mat B;
71: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
72: Mat_SeqBAIJ *b;
73: PetscInt *ai=a->i,m=A->rmap->N,n=A->cmap->N,i,*rowlengths,bs=PetscAbs(A->rmap->bs);
75: if (reuse != MAT_REUSE_MATRIX) {
76: PetscMalloc1(m/bs,&rowlengths);
77: for (i=0; i<m/bs; i++) {
78: rowlengths[i] = (ai[i*bs+1] - ai[i*bs])/bs;
79: }
80: MatCreate(PetscObjectComm((PetscObject)A),&B);
81: MatSetSizes(B,m,n,m,n);
82: MatSetType(B,MATSEQBAIJ);
83: MatSeqBAIJSetPreallocation(B,bs,0,rowlengths);
84: PetscFree(rowlengths);
85: } else B = *newmat;
87: if (bs == 1) {
88: b = (Mat_SeqBAIJ*)(B->data);
90: PetscArraycpy(b->i,a->i,m+1);
91: PetscArraycpy(b->ilen,a->ilen,m);
92: PetscArraycpy(b->j,a->j,a->nz);
93: PetscArraycpy(b->a,a->a,a->nz);
95: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
96: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
98: } else {
99: /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
100: /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
101: /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
102: MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);
103: }
105: if (reuse == MAT_INPLACE_MATRIX) {
106: MatHeaderReplace(A,&B);
107: } else *newmat = B;
108: return 0;
109: }