Actual source code: mcrl.c
2: /*
3: Defines a matrix-vector product for the MATMPIAIJCRL matrix class.
4: This class is derived from the MATMPIAIJ class and retains the
5: compressed row storage (aka Yale sparse matrix format) but augments
6: it with a column oriented storage that is more efficient for
7: matrix vector products on Vector machines.
9: CRL stands for constant row length (that is the same number of columns
10: is kept (padded with zeros) for each row of the sparse matrix.
12: See src/mat/impls/aij/seq/crl/crl.c for the sequential version
13: */
15: #include <../src/mat/impls/aij/mpi/mpiaij.h>
16: #include <../src/mat/impls/aij/seq/crl/crl.h>
18: PetscErrorCode MatDestroy_MPIAIJCRL(Mat A)
19: {
20: Mat_AIJCRL *aijcrl = (Mat_AIJCRL*) A->spptr;
22: if (aijcrl) {
23: PetscFree2(aijcrl->acols,aijcrl->icols);
24: VecDestroy(&aijcrl->fwork);
25: VecDestroy(&aijcrl->xwork);
26: PetscFree(aijcrl->array);
27: }
28: PetscFree(A->spptr);
30: PetscObjectChangeTypeName((PetscObject)A, MATMPIAIJ);
31: MatDestroy_MPIAIJ(A);
32: return 0;
33: }
35: PetscErrorCode MatMPIAIJCRL_create_aijcrl(Mat A)
36: {
37: Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A)->data;
38: Mat_SeqAIJ *Aij = (Mat_SeqAIJ*)(a->A->data), *Bij = (Mat_SeqAIJ*)(a->B->data);
39: Mat_AIJCRL *aijcrl = (Mat_AIJCRL*) A->spptr;
40: PetscInt m = A->rmap->n; /* Number of rows in the matrix. */
41: PetscInt nd = a->A->cmap->n; /* number of columns in diagonal portion */
42: PetscInt *aj = Aij->j,*bj = Bij->j; /* From the CSR representation; points to the beginning of each row. */
43: PetscInt i, j,rmax = 0,*icols, *ailen = Aij->ilen, *bilen = Bij->ilen;
44: PetscScalar *aa = Aij->a,*ba = Bij->a,*acols,*array;
46: /* determine the row with the most columns */
47: for (i=0; i<m; i++) {
48: rmax = PetscMax(rmax,ailen[i]+bilen[i]);
49: }
50: aijcrl->nz = Aij->nz+Bij->nz;
51: aijcrl->m = A->rmap->n;
52: aijcrl->rmax = rmax;
54: PetscFree2(aijcrl->acols,aijcrl->icols);
55: PetscMalloc2(rmax*m,&aijcrl->acols,rmax*m,&aijcrl->icols);
56: acols = aijcrl->acols;
57: icols = aijcrl->icols;
58: for (i=0; i<m; i++) {
59: for (j=0; j<ailen[i]; j++) {
60: acols[j*m+i] = *aa++;
61: icols[j*m+i] = *aj++;
62: }
63: for (; j<ailen[i]+bilen[i]; j++) {
64: acols[j*m+i] = *ba++;
65: icols[j*m+i] = nd + *bj++;
66: }
67: for (; j<rmax; j++) { /* empty column entries */
68: acols[j*m+i] = 0.0;
69: icols[j*m+i] = (j) ? icols[(j-1)*m+i] : 0; /* handle case where row is EMPTY */
70: }
71: }
72: PetscInfo(A,"Percentage of 0's introduced for vectorized multiply %g\n",1.0-((double)(aijcrl->nz))/((double)(rmax*m)));
74: PetscFree(aijcrl->array);
75: PetscMalloc1(a->B->cmap->n+nd,&array);
76: /* xwork array is actually B->n+nd long, but we define xwork this length so can copy into it */
77: VecDestroy(&aijcrl->xwork);
78: VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,nd,PETSC_DECIDE,array,&aijcrl->xwork);
79: VecDestroy(&aijcrl->fwork);
80: VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->B->cmap->n,array+nd,&aijcrl->fwork);
82: aijcrl->array = array;
83: aijcrl->xscat = a->Mvctx;
84: return 0;
85: }
87: PetscErrorCode MatAssemblyEnd_MPIAIJCRL(Mat A, MatAssemblyType mode)
88: {
89: Mat_MPIAIJ *a = (Mat_MPIAIJ*)A->data;
90: Mat_SeqAIJ *Aij = (Mat_SeqAIJ*)(a->A->data), *Bij = (Mat_SeqAIJ*)(a->A->data);
92: Aij->inode.use = PETSC_FALSE;
93: Bij->inode.use = PETSC_FALSE;
95: MatAssemblyEnd_MPIAIJ(A,mode);
96: if (mode == MAT_FLUSH_ASSEMBLY) return 0;
98: /* Now calculate the permutation and grouping information. */
99: MatMPIAIJCRL_create_aijcrl(A);
100: return 0;
101: }
103: extern PetscErrorCode MatMult_AIJCRL(Mat,Vec,Vec);
104: extern PetscErrorCode MatDuplicate_AIJCRL(Mat,MatDuplicateOption,Mat*);
106: /* MatConvert_MPIAIJ_MPIAIJCRL converts a MPIAIJ matrix into a
107: * MPIAIJCRL matrix. This routine is called by the MatCreate_MPIAIJCRL()
108: * routine, but can also be used to convert an assembled MPIAIJ matrix
109: * into a MPIAIJCRL one. */
111: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_MPIAIJCRL(Mat A,MatType type,MatReuse reuse,Mat *newmat)
112: {
113: Mat B = *newmat;
114: Mat_AIJCRL *aijcrl;
116: if (reuse == MAT_INITIAL_MATRIX) {
117: MatDuplicate(A,MAT_COPY_VALUES,&B);
118: }
120: PetscNewLog(B,&aijcrl);
121: B->spptr = (void*) aijcrl;
123: /* Set function pointers for methods that we inherit from AIJ but override. */
124: B->ops->duplicate = MatDuplicate_AIJCRL;
125: B->ops->assemblyend = MatAssemblyEnd_MPIAIJCRL;
126: B->ops->destroy = MatDestroy_MPIAIJCRL;
127: B->ops->mult = MatMult_AIJCRL;
129: /* If A has already been assembled, compute the permutation. */
130: if (A->assembled) {
131: MatMPIAIJCRL_create_aijcrl(B);
132: }
133: PetscObjectChangeTypeName((PetscObject)B,MATMPIAIJCRL);
134: *newmat = B;
135: return 0;
136: }
138: /*@C
139: MatCreateMPIAIJCRL - Creates a sparse matrix of type MPIAIJCRL.
140: This type inherits from AIJ, but stores some additional
141: information that is used to allow better vectorization of
142: the matrix-vector product. At the cost of increased storage, the AIJ formatted
143: matrix can be copied to a format in which pieces of the matrix are
144: stored in ELLPACK format, allowing the vectorized matrix multiply
145: routine to use stride-1 memory accesses. As with the AIJ type, it is
146: important to preallocate matrix storage in order to get good assembly
147: performance.
149: Collective
151: Input Parameters:
152: + comm - MPI communicator, set to PETSC_COMM_SELF
153: . m - number of rows
154: . n - number of columns
155: . nz - number of nonzeros per row (same for all rows)
156: - nnz - array containing the number of nonzeros in the various rows
157: (possibly different for each row) or NULL
159: Output Parameter:
160: . A - the matrix
162: Notes:
163: If nnz is given then nz is ignored
165: Level: intermediate
167: .seealso: MatCreate(), MatCreateMPIAIJPERM(), MatSetValues()
168: @*/
169: PetscErrorCode MatCreateMPIAIJCRL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],PetscInt onz,const PetscInt onnz[],Mat *A)
170: {
171: MatCreate(comm,A);
172: MatSetSizes(*A,m,n,m,n);
173: MatSetType(*A,MATMPIAIJCRL);
174: MatMPIAIJSetPreallocation_MPIAIJ(*A,nz,(PetscInt*)nnz,onz,(PetscInt*)onnz);
175: return 0;
176: }
178: PETSC_EXTERN PetscErrorCode MatCreate_MPIAIJCRL(Mat A)
179: {
180: MatSetType(A,MATMPIAIJ);
181: MatConvert_MPIAIJ_MPIAIJCRL(A,MATMPIAIJCRL,MAT_INPLACE_MATRIX,&A);
182: return 0;
183: }