Actual source code: convert.c
petsc-3.5.4 2015-05-23
2: #include <petsc-private/matimpl.h>
6: /*
7: MatConvert_Basic - Converts from any input format to another format. For
8: parallel formats, the new matrix distribution is determined by PETSc.
10: Does not do preallocation so in general will be slow
11: */
12: PetscErrorCode MatConvert_Basic(Mat mat, MatType newtype,MatReuse reuse,Mat *newmat)
13: {
14: Mat M;
15: const PetscScalar *vwork;
16: PetscErrorCode ierr;
17: PetscInt i,j,nz,m,n,rstart,rend,lm,ln,prbs,pcbs,cstart,cend,*dnz,*onz;
18: const PetscInt *cwork;
19: PetscBool isseqsbaij,ismpisbaij,isseqbaij,ismpibaij,isseqdense,ismpidense;
22: MatGetSize(mat,&m,&n);
23: MatGetLocalSize(mat,&lm,&ln);
25: if (ln == n) ln = PETSC_DECIDE; /* try to preserve column ownership */
27: MatCreate(PetscObjectComm((PetscObject)mat),&M);
28: MatSetSizes(M,lm,ln,m,n);
29: MatSetBlockSizesFromMats(M,mat,mat);
30: MatSetType(M,newtype);
31: MatGetOwnershipRange(mat,&rstart,&rend);
33: PetscObjectTypeCompare((PetscObject)M,MATSEQSBAIJ,&isseqsbaij);
34: PetscObjectTypeCompare((PetscObject)M,MATMPISBAIJ,&ismpisbaij);
35: if (isseqsbaij || ismpisbaij) {MatSetOption(M,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);}
36: PetscObjectTypeCompare((PetscObject)M,MATSEQBAIJ,&isseqbaij);
37: PetscObjectTypeCompare((PetscObject)M,MATMPIBAIJ,&ismpibaij);
38: PetscObjectTypeCompare((PetscObject)M,MATSEQDENSE,&isseqdense);
39: PetscObjectTypeCompare((PetscObject)M,MATMPIDENSE,&ismpidense);
41: if (isseqdense) {
42: MatSeqDenseSetPreallocation(M,NULL);
43: } else if (ismpidense) {
44: MatMPIDenseSetPreallocation(M,NULL);
45: } else {
46: /* Preallocation block sizes. (S)BAIJ matrices will have one index per block. */
47: prbs = (isseqbaij || ismpibaij || isseqsbaij || ismpisbaij) ? PetscAbs(M->rmap->bs) : 1;
48: pcbs = (isseqbaij || ismpibaij || isseqsbaij || ismpisbaij) ? PetscAbs(M->cmap->bs) : 1;
50: PetscMalloc2(lm/prbs,&dnz,lm/prbs,&onz);
51: MatGetOwnershipRangeColumn(mat,&cstart,&cend);
52: for (i=0; i<lm; i+=prbs) {
53: MatGetRow(mat,rstart+i,&nz,&cwork,NULL);
54: dnz[i] = 0;
55: onz[i] = 0;
56: for (j=0; j<nz; j+=pcbs) {
57: if ((isseqsbaij || ismpisbaij) && cwork[j] < rstart+i) continue;
58: if (cstart <= cwork[j] && cwork[j] < cend) dnz[i]++;
59: else onz[i]++;
60: }
61: MatRestoreRow(mat,rstart+i,&nz,&cwork,NULL);
62: }
63: MatXAIJSetPreallocation(M,PETSC_DECIDE,dnz,onz,dnz,onz);
64: PetscFree2(dnz,onz);
65: }
67: for (i=rstart; i<rend; i++) {
68: MatGetRow(mat,i,&nz,&cwork,&vwork);
69: MatSetValues(M,1,&i,nz,cwork,vwork,INSERT_VALUES);
70: MatRestoreRow(mat,i,&nz,&cwork,&vwork);
71: }
72: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
75: if (reuse == MAT_REUSE_MATRIX) {
76: MatHeaderReplace(mat,M);
77: } else {
78: *newmat = M;
79: }
80: return(0);
81: }