Actual source code: aijbaij.c
petsc-3.8.4 2018-03-24
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;
11: PetscInt bs = A->rmap->bs,*ai = a->i,*aj = a->j,n = A->rmap->N/bs,i,j,k;
12: PetscInt *rowlengths,*rows,*cols,maxlen = 0,ncols;
13: MatScalar *aa = a->a;
16: if (reuse == MAT_REUSE_MATRIX) {
17: B = *newmat;
18: for (i=0; i<n; i++) {
19: maxlen = PetscMax(maxlen,(ai[i+1] - ai[i]));
20: }
21: } else {
22: PetscMalloc1(n*bs,&rowlengths);
23: for (i=0; i<n; i++) {
24: maxlen = PetscMax(maxlen,(ai[i+1] - ai[i]));
25: for (j=0; j<bs; j++) {
26: rowlengths[i*bs+j] = bs*(ai[i+1] - ai[i]);
27: }
28: }
29: MatCreate(PetscObjectComm((PetscObject)A),&B);
30: MatSetType(B,MATSEQAIJ);
31: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
32: MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
33: MatSeqAIJSetPreallocation(B,0,rowlengths);
34: PetscFree(rowlengths);
35: }
36: b = (Mat_SeqAIJ*)B->data;
37: roworiented = b->roworiented;
39: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);
40: PetscMalloc1(bs,&rows);
41: PetscMalloc1(bs*maxlen,&cols);
42: for (i=0; i<n; i++) {
43: for (j=0; j<bs; j++) {
44: rows[j] = i*bs+j;
45: }
46: ncols = ai[i+1] - ai[i];
47: for (k=0; k<ncols; k++) {
48: for (j=0; j<bs; j++) {
49: cols[k*bs+j] = bs*(*aj) + j;
50: }
51: aj++;
52: }
53: MatSetValues(B,bs,rows,bs*ncols,cols,aa,INSERT_VALUES);
54: aa += ncols*bs*bs;
55: }
56: PetscFree(cols);
57: PetscFree(rows);
58: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
59: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
60: MatSetOption(B,MAT_ROW_ORIENTED,roworiented);
62: if (reuse == MAT_INPLACE_MATRIX) {
63: MatHeaderReplace(A,&B);
64: } else {
65: *newmat = B;
66: }
67: return(0);
68: }
70: #include <../src/mat/impls/aij/seq/aij.h>
72: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
73: {
74: Mat B;
75: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
76: Mat_SeqBAIJ *b;
78: PetscInt *ai=a->i,m=A->rmap->N,n=A->cmap->N,i,*rowlengths;
81: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
82: if (A->rmap->bs > 1) {
83: MatConvert_Basic(A,newtype,reuse,newmat);
84: return(0);
85: }
87: PetscMalloc1(m,&rowlengths);
88: for (i=0; i<m; i++) {
89: rowlengths[i] = ai[i+1] - ai[i];
90: }
91: MatCreate(PetscObjectComm((PetscObject)A),&B);
92: MatSetSizes(B,m,n,m,n);
93: MatSetType(B,MATSEQBAIJ);
94: MatSeqBAIJSetPreallocation(B,1,0,rowlengths);
95: PetscFree(rowlengths);
97: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
99: b = (Mat_SeqBAIJ*)(B->data);
101: PetscMemcpy(b->i,a->i,(m+1)*sizeof(PetscInt));
102: PetscMemcpy(b->ilen,a->ilen,m*sizeof(PetscInt));
103: PetscMemcpy(b->j,a->j,a->nz*sizeof(PetscInt));
104: PetscMemcpy(b->a,a->a,a->nz*sizeof(MatScalar));
106: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
107: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
109: if (reuse == MAT_INPLACE_MATRIX) {
110: MatHeaderReplace(A,&B);
111: } else {
112: *newmat = B;
113: }
114: return(0);
115: }