Actual source code: mpisbstream.c
petsc-3.6.1 2015-08-06
1: #define PETSCMAT_DLL
3: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
4: #include <../src/mat/impls/sbaij/seq/sbstream/sbstream.h>
6: extern PetscErrorCode MatMult_SeqSBSTRM_4(Mat,Vec,Vec);
7: extern PetscErrorCode MatMult_SeqSBSTRM_5(Mat,Vec,Vec);
8: extern PetscErrorCode MatMultTranspose_SeqBSTRM_4(Mat,Vec,Vec);
9: extern PetscErrorCode MatMultTranspose_SeqBSTRM_5(Mat,Vec,Vec);
10: extern PetscErrorCode MatMultAdd_SeqSBSTRM_4(Mat,Vec,Vec,Vec);
11: extern PetscErrorCode MatMultAdd_SeqSBSTRM_5(Mat,Vec,Vec,Vec);
12: extern PetscErrorCode MatMultAdd_SeqBSTRM_4(Mat,Vec,Vec,Vec);
13: extern PetscErrorCode MatMultAdd_SeqBSTRM_5(Mat,Vec,Vec,Vec);
17: PetscErrorCode MPISBSTRM_create_sbstrm(Mat A)
18: {
19: Mat_MPISBAIJ *a = (Mat_MPISBAIJ*)A->data;
20: Mat_SeqSBAIJ *Aij = (Mat_SeqSBAIJ*)(a->A->data), *Bij = (Mat_SeqSBAIJ*)(a->B->data);
21: /*
22: */
23: Mat_SeqSBSTRM *sbstrmA, *sbstrmB;
24: PetscInt MROW = Aij->mbs, bs = a->A->rmap->bs;
26: /* PetscInt m = A->rmap->n;*/ /* Number of rows in the matrix. */
27: /* PetscInt nd = a->A->cmap->n;*/ /* number of columns in diagonal portion */
28: PetscInt *ai = Aij->i, *bi = Bij->i; /* From the CSR representation; points to the beginning of each row. */
29: PetscInt i,j,k;
30: PetscScalar *aa = Aij->a,*ba = Bij->a;
32: PetscInt bs2, rbs, cbs, slen, blen;
34: PetscScalar **asp;
35: PetscScalar **bsp;
38: /* printf(" --- in MPISBSTRM_create_sbstrm, m=%d, nd=%d, bs=%d, MROW=%d\n", m,nd,bs,MROW); */
40: rbs = cbs = bs;
41: bs2 = bs*bs;
42: blen = ai[MROW]-ai[0];
43: slen = blen*bs;
45: /* printf(" --- blen=%d, slen=%d\n", blen, slen); */
47: PetscNewLog(a->A,&sbstrmA);
48: a->A->spptr = (void*) sbstrmA;
49: sbstrmA = (Mat_SeqSBSTRM*) a->A->spptr;
50: sbstrmA->rbs = sbstrmA->cbs = bs;
51: PetscMalloc1(bs2*blen, &sbstrmA->as);
53: PetscMalloc1(rbs, &asp);
55: for (i=0; i<rbs; i++) asp[i] = sbstrmA->as + i*slen;
57: for (k=0; k<blen; k++) {
58: for (j=0; j<cbs; j++) {
59: for (i=0; i<rbs; i++) asp[i][k*cbs+j] = aa[k*bs2+j*rbs+i];
60: }
61: }
62: switch (bs) {
63: case 4:
64: a->A->ops->mult = MatMult_SeqSBSTRM_4;
65: a->A->ops->multtranspose = MatMult_SeqSBSTRM_4;
66: /* a->A->ops->sor = MatSOR_SeqSBSTRM_4; */
67: break;
68: case 5:
69: a->A->ops->mult = MatMult_SeqSBSTRM_5;
70: a->A->ops->multtranspose = MatMult_SeqSBSTRM_5;
71: /* a->A->ops->sor = MatSOR_SeqSBSTRM_5; */
72: break;
73: default:
74: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"not supported for block size %D yet",bs);
75: }
76: PetscFree(asp);
79: /*.....*/
80: blen = bi[MROW]-bi[0];
81: slen = blen*bs;
82: PetscNewLog(a->B,&sbstrmB);
83: a->B->spptr = (void*) sbstrmB;
84: sbstrmB = (Mat_SeqSBSTRM*) a->B->spptr;
85: sbstrmB->rbs = sbstrmB->cbs = bs;
87: PetscMalloc1(bs2*blen, &sbstrmB->as);
89: PetscMalloc1(rbs, &bsp);
91: for (i=0; i<rbs; i++) bsp[i] = sbstrmB->as + i*slen;
93: for (k=0; k<blen; k++) {
94: for (j=0; j<cbs; j++) {
95: for (i=0; i<rbs; i++) bsp[i][k*cbs+j] = ba[k*bs2+j*rbs+i];
96: }
97: }
98: switch (bs) {
99: case 4:
100: /* a->B->ops->mult = MatMult_SeqSBSTRM_4; */
101: a->B->ops->multtranspose = MatMultTranspose_SeqBSTRM_4;
102: a->B->ops->multadd = MatMultAdd_SeqBSTRM_4;
103: /* a->B->ops->multtransposeadd = MatMultAdd_SeqSBSTRM_4; */
104: break;
105: case 5:
106: /* a->B->ops->mult = MatMult_SeqSBSTRM_5; */
107: a->B->ops->multtranspose = MatMultTranspose_SeqBSTRM_5;
108: a->B->ops->multadd = MatMultAdd_SeqBSTRM_5;
109: /* a->B->ops->multtransposeadd = MatMultAdd_SeqSBSTRM_5; */
110: break;
111: default:
112: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"not supported for block size %D yet",bs);
113: }
114: PetscFree(bsp);
115: return(0);
116: }
119: extern PetscErrorCode MatAssemblyEnd_MPISBAIJ(Mat,MatAssemblyType);
122: PetscErrorCode MatAssemblyEnd_MPISBSTRM(Mat A, MatAssemblyType mode)
123: {
127: /*
128: Aij->inode.use = PETSC_FALSE;
129: Bij->inode.use = PETSC_FALSE;
130: */
131: MatAssemblyEnd_MPISBAIJ(A,mode);
132: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
134: /* Now calculate the permutation and grouping information. */
135: MPISBSTRM_create_sbstrm(A);
136: return(0);
137: }
142: PetscErrorCode MatCreateMPISBSTRM(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[],Mat *A)
143: {
145: PetscMPIInt size;
148: MatCreate(comm,A);
149: MatSetSizes(*A,m,n,M,N);
150: MPI_Comm_size(comm,&size);
151: if (size > 1) {
152: MatSetType(*A,MATMPISBSTRM);
153: MatMPISBAIJSetPreallocation(*A,bs,d_nz,d_nnz,o_nz,o_nnz);
154: } else {
155: MatSetType(*A,MATSEQSBSTRM);
156: MatSeqSBAIJSetPreallocation(*A,bs,d_nz,d_nnz);
157: }
158: return(0);
159: }
161: PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat,MatType,MatReuse,Mat*);
162: extern PetscErrorCode MatMPISBAIJSetPreallocation_MPISBAIJ(Mat,PetscInt,PetscInt,const PetscInt *,PetscInt,const PetscInt*);
166: PetscErrorCode MatMPISBAIJSetPreallocation_MPISBSTRM(Mat B,PetscInt bs,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
167: {
171: MatMPISBAIJSetPreallocation_MPISBAIJ(B,bs,d_nz,d_nnz,o_nz,o_nnz);
172: return(0);
173: }
177: PETSC_EXTERN PetscErrorCode MatConvert_MPISBAIJ_MPISBSTRM(Mat A,MatType type,MatReuse reuse,Mat *newmat)
178: {
180: Mat B = *newmat;
181: Mat_SeqSBSTRM *sbstrm;
184: if (reuse == MAT_INITIAL_MATRIX) {
185: MatDuplicate(A,MAT_COPY_VALUES,&B);
186: }
187: /* printf(" --- in MatConvert_MPISBAIJ_MPISBSTRM -- 1 \n"); */
189: PetscNewLog(B,&sbstrm);
190: B->spptr = (void*)sbstrm;
192: /* Set function pointers for methods that we inherit from AIJ but override.
193: B->ops->duplicate = MatDuplicate_SBSTRM;
194: B->ops->mult = MatMult_SBSTRM;
195: B->ops->destroy = MatDestroy_MPISBSTRM;
196: */
197: B->ops->assemblyend = MatAssemblyEnd_MPISBSTRM;
199: /* If A has already been assembled, compute the permutation. */
200: if (A->assembled) {
201: MPISBSTRM_create_sbstrm(B);
202: }
204: PetscObjectChangeTypeName((PetscObject) B, MATMPISBSTRM);
205: PetscObjectComposeFunction((PetscObject)B,"MatMPISBAIJSetPreallocation_C",MatMPISBAIJSetPreallocation_MPISBSTRM);
206: *newmat = B;
207: return(0);
208: }
212: PETSC_EXTERN PetscErrorCode MatCreate_MPISBSTRM(Mat A)
213: {
217: MatSetType(A,MATMPISBAIJ);
218: MatConvert_MPISBAIJ_MPISBSTRM(A,MATMPISBSTRM,MAT_REUSE_MATRIX,&A);
219: return(0);
220: }
224: PETSC_EXTERN PetscErrorCode MatCreate_SBSTRM(Mat A)
225: {
227: PetscMPIInt size;
230: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
231: if (size == 1) {
232: MatSetType(A,MATSEQSBSTRM);
233: } else {
234: MatSetType(A,MATMPISBSTRM);
235: }
236: return(0);
237: }