Actual source code: aijbaij.c
petsc-3.7.3 2016-08-01
2: #include <../src/mat/impls/baij/seq/baij.h>
6: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
7: {
8: Mat B;
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: PetscMalloc1(n*bs,&rowlengths);
17: for (i=0; i<n; i++) {
18: maxlen = PetscMax(maxlen,(ai[i+1] - ai[i]));
19: for (j=0; j<bs; j++) {
20: rowlengths[i*bs+j] = bs*(ai[i+1] - ai[i]);
21: }
22: }
23: MatCreate(PetscObjectComm((PetscObject)A),&B);
24: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
25: MatSetType(B,MATSEQAIJ);
26: MatSeqAIJSetPreallocation(B,0,rowlengths);
27: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);
28: PetscFree(rowlengths);
30: PetscMalloc1(bs,&rows);
31: PetscMalloc1(bs*maxlen,&cols);
32: for (i=0; i<n; i++) {
33: for (j=0; j<bs; j++) {
34: rows[j] = i*bs+j;
35: }
36: ncols = ai[i+1] - ai[i];
37: for (k=0; k<ncols; k++) {
38: for (j=0; j<bs; j++) {
39: cols[k*bs+j] = bs*(*aj) + j;
40: }
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);
51: B->rmap->bs = A->rmap->bs;
53: if (reuse == MAT_INPLACE_MATRIX) {
54: MatHeaderReplace(A,&B);
55: } else {
56: *newmat = B;
57: }
58: return(0);
59: }
61: #include <../src/mat/impls/aij/seq/aij.h>
65: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
66: {
67: Mat B;
68: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
69: Mat_SeqBAIJ *b;
71: PetscInt *ai=a->i,m=A->rmap->N,n=A->cmap->N,i,*rowlengths;
74: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
75: if (A->rmap->bs > 1) {
76: MatConvert_Basic(A,newtype,reuse,newmat);
77: return(0);
78: }
80: PetscMalloc1(m,&rowlengths);
81: for (i=0; i<m; i++) {
82: rowlengths[i] = ai[i+1] - ai[i];
83: }
84: MatCreate(PetscObjectComm((PetscObject)A),&B);
85: MatSetSizes(B,m,n,m,n);
86: MatSetType(B,MATSEQBAIJ);
87: MatSeqBAIJSetPreallocation_SeqBAIJ(B,1,0,rowlengths);
88: PetscFree(rowlengths);
90: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
92: b = (Mat_SeqBAIJ*)(B->data);
94: PetscMemcpy(b->i,a->i,(m+1)*sizeof(PetscInt));
95: PetscMemcpy(b->ilen,a->ilen,m*sizeof(PetscInt));
96: PetscMemcpy(b->j,a->j,a->nz*sizeof(PetscInt));
97: PetscMemcpy(b->a,a->a,a->nz*sizeof(MatScalar));
99: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
100: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
102: if (reuse == MAT_INPLACE_MATRIX) {
103: MatHeaderReplace(A,&B);
104: } else {
105: *newmat = B;
106: }
107: return(0);
108: }