1 #ifndef lint 2 static char vcid[] = "$Id: convert.c,v 1.42 1996/01/12 22:31:53 bsmith Exp bsmith $"; 3 #endif 4 5 #include "mpiaij.h" 6 #include "mpibdiag.h" 7 8 /* This file contains a generic conversion routine and implementation specific 9 versions for increased efficiency. */ 10 11 /* 12 MatConvert_Basic - Converts from any input format to another format. For 13 parallel formats, the new matrix distribution is determined by PETSc. 14 */ 15 int MatConvert_Basic(Mat mat,MatType newtype,Mat *M) 16 { 17 Scalar *vwork; 18 int ierr, i, nz, m, n, *cwork, rstart, rend,flg; 19 20 ierr = MatGetSize(mat,&m,&n); CHKERRQ(ierr); 21 if (newtype == MATSAME) newtype = (MatType)mat->type; 22 switch (newtype) { 23 case MATSEQAIJ: 24 ierr = MatCreateSeqAIJ(mat->comm,m,n,0,PETSC_NULL,M); CHKERRQ(ierr); 25 break; 26 case MATMPIROWBS: 27 if (m != n) SETERRQ(1,"MatConvert:MATMPIROWBS matrix must be square"); 28 ierr = MatCreateMPIRowbs(MPI_COMM_WORLD,PETSC_DECIDE,m,0,PETSC_NULL, 29 PETSC_NULL,M); CHKERRQ(ierr); 30 break; 31 case MATMPIAIJ: 32 ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 33 m,n,0,PETSC_NULL,0,PETSC_NULL,M); CHKERRQ(ierr); 34 break; 35 case MATSEQDENSE: 36 ierr = MatCreateSeqDense(mat->comm,m,n,PETSC_NULL,M); CHKERRQ(ierr); 37 break; 38 case MATMPIDENSE: 39 ierr = MatCreateMPIDense(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 40 m,n,PETSC_NULL,M); CHKERRQ(ierr); 41 break; 42 case MATSEQBDIAG: 43 { 44 int nb = 1; /* Default block size = 1 */ 45 ierr = OptionsGetInt(PETSC_NULL,"-mat_bdiag_bsize",&nb,&flg); CHKERRQ(ierr); 46 ierr = MatCreateSeqBDiag(mat->comm,m,n,0,nb,PETSC_NULL,PETSC_NULL,M); CHKERRQ(ierr); 47 break; 48 } 49 case MATMPIBDIAG: 50 { 51 int nb = 1; /* Default block size = 1 */ 52 ierr = OptionsGetInt(PETSC_NULL,"-mat_bdiag_bsize",&nb,&flg); CHKERRQ(ierr); 53 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,0,nb,PETSC_NULL, 54 PETSC_NULL,M); CHKERRQ(ierr); 55 break; 56 } 57 default: 58 SETERRQ(1,"MatConvert:Matrix type is not currently supported"); 59 } 60 ierr = MatGetOwnershipRange(*M,&rstart,&rend); CHKERRQ(ierr); 61 for (i=rstart; i<rend; i++) { 62 ierr = MatGetRow(mat,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 63 ierr = MatSetValues(*M,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 64 ierr = MatRestoreRow(mat,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 65 } 66 ierr = MatAssemblyBegin(*M,FINAL_ASSEMBLY); CHKERRQ(ierr); 67 ierr = MatAssemblyEnd(*M,FINAL_ASSEMBLY); CHKERRQ(ierr); 68 return 0; 69 } 70 /* -------------------------------------------------------------- */ 71 /* 72 MatConvert_SeqAIJ - Converts from MATSEQAIJ format to another format. For 73 parallel formats, the new matrix distribution is determined by PETSc. 74 */ 75 int MatConvert_SeqAIJ(Mat A, MatType newtype, Mat *B) 76 { 77 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 78 Scalar *vwork; 79 int i, ierr, nz, m = a->m, n = a->n, *cwork, rstart, rend,flg; 80 81 switch (newtype) { 82 case MATMPIROWBS: 83 if (m != n) SETERRQ(1,"MatConvert_SeqAIJ:MATMPIROWBS matrix must be square"); 84 ierr = MatCreateMPIRowbs(MPI_COMM_WORLD,PETSC_DECIDE,m,0,PETSC_NULL,PETSC_NULL,B); 85 CHKERRQ(ierr); 86 break; 87 case MATMPIAIJ: 88 ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 89 m,n,0,PETSC_NULL,0,PETSC_NULL,B); CHKERRQ(ierr); 90 break; 91 case MATSEQDENSE: 92 ierr = MatCreateSeqDense(A->comm,m,n,PETSC_NULL,B); CHKERRQ(ierr); 93 break; 94 case MATMPIDENSE: 95 ierr = MatCreateMPIDense(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 96 m,n,PETSC_NULL,B); CHKERRQ(ierr); 97 break; 98 case MATSEQBDIAG: 99 { 100 int nb = 1; /* Default block size = 1 */ 101 ierr = OptionsGetInt(PETSC_NULL,"-mat_bdiag_bsize",&nb,&flg); CHKERRQ(ierr); 102 ierr = MatCreateSeqBDiag(A->comm,m,n,0,nb,PETSC_NULL,PETSC_NULL,B); CHKERRQ(ierr); 103 break; 104 } 105 case MATMPIBDIAG: 106 { 107 int nb = 1; /* Default block size = 1 */ 108 ierr = OptionsGetInt(PETSC_NULL,"-mat_bdiag_bsize",&nb,&flg); CHKERRQ(ierr); 109 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,0,nb,PETSC_NULL, 110 PETSC_NULL,B); CHKERRQ(ierr); 111 break; 112 } 113 default: 114 SETERRQ(1,"MatConvert_SeqAIJ:Matrix type is not currently supported"); 115 } 116 ierr = MatGetOwnershipRange(*B,&rstart,&rend); CHKERRQ(ierr); 117 for (i=rstart; i<rend; i++) { 118 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 119 ierr = MatSetValues(*B,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 120 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 121 } 122 ierr = MatAssemblyBegin(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 123 ierr = MatAssemblyEnd(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 124 return 0; 125 } 126 /* ------------------------------------------------------------------ */ 127 /* 128 MatConvert_MPIAIJ - Converts from MATMPIAIJ format to another 129 parallel format. 130 */ 131 int MatConvert_MPIAIJ(Mat A, MatType newtype, Mat *B) 132 { 133 Mat_MPIAIJ *a = (Mat_MPIAIJ *) A->data; 134 int ierr, nz, i, ig, rstart = a->rstart, m = a->m, *cwork; 135 Scalar *vwork; 136 137 SETERRQ(1,"MatConvert_MPIAIJ:Not currently suported"); 138 139 /* Each processor converts its local rows */ 140 for (i=0; i<m; i++) { 141 ig = i + rstart; 142 ierr = MatGetRow(A,ig,&nz,&cwork,&vwork); CHKERRQ(ierr); 143 ierr = MatSetValues(*B,1,&ig,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 144 ierr = MatRestoreRow(A,ig,&nz,&cwork,&vwork); CHKERRQ(ierr); 145 } 146 ierr = MatAssemblyBegin(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 147 ierr = MatAssemblyEnd(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 148 return 0; 149 } 150 /* ------------------------------------------------------------------ */ 151 /* 152 MatConvert_SeqBDiag - Converts from MATSEQBDiag format to another format. For 153 parallel formats, the new matrix distribution is determined by PETSc. 154 */ 155 int MatConvert_SeqBDiag(Mat A, MatType newtype, Mat *B) 156 { 157 Mat_SeqBDiag *a = (Mat_SeqBDiag *) A->data; 158 Scalar *vwork, *vw2; 159 int i, ierr, nz, m = a->m, n = a->n, *cwork, rstart, rend; 160 int j, *cw2, ict; 161 162 /* rough over-estimate; could refine for individual rows */ 163 nz = PetscMin(n,a->nd*a->nb); 164 switch (newtype) { 165 case MATSEQAIJ: 166 ierr = MatCreateSeqAIJ(A->comm,m,n,nz,PETSC_NULL,B); CHKERRQ(ierr); 167 break; 168 case MATMPIROWBS: 169 if (m != n) SETERRQ(1,"MatConvert_SeqBDiag:MATMPIROWBS matrix must be square"); 170 ierr = MatCreateMPIRowbs(MPI_COMM_WORLD,PETSC_DECIDE,m,0,PETSC_NULL,PETSC_NULL, 171 B); CHKERRQ(ierr); 172 break; 173 case MATMPIAIJ: 174 ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 175 m,n,0,PETSC_NULL,0,PETSC_NULL,B); CHKERRQ(ierr); 176 break; 177 case MATSEQDENSE: 178 ierr = MatCreateSeqDense(A->comm,m,n,PETSC_NULL,B); CHKERRQ(ierr); 179 break; 180 case MATMPIDENSE: 181 ierr = MatCreateMPIDense(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 182 m,n,PETSC_NULL,B); CHKERRQ(ierr); 183 break; 184 case MATMPIBDIAG: 185 { 186 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,a->nd,a->nb, 187 PETSC_NULL,PETSC_NULL,B); CHKERRQ(ierr); 188 break; 189 } 190 default: 191 SETERRQ(1,"MatConvert_SeqBDiag:Matrix type is not currently supported"); 192 } 193 ierr = MatGetOwnershipRange(*B,&rstart,&rend); CHKERRQ(ierr); 194 195 cw2 = (int *)PetscMalloc( n * sizeof(int) ); CHKPTRQ(cw2); 196 vw2 = (Scalar *)PetscMalloc( n * sizeof(Scalar) ); CHKPTRQ(vw2); 197 for (i=rstart; i<rend; i++) { 198 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 199 ict = 0; /* strip out the zero elements ... is this what we really want? */ 200 for (j=0; j<nz; j++) { 201 if (vwork[j] != 0) {vw2[ict] = vwork[j]; cw2[ict] = cwork[j]; ict++;} 202 } 203 if (ict) 204 {ierr = MatSetValues(*B,1,&i,ict,cw2,vw2,INSERT_VALUES); CHKERRQ(ierr);} 205 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 206 } 207 PetscFree(cw2); PetscFree(vw2); 208 ierr = MatAssemblyBegin(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 209 ierr = MatAssemblyEnd(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 210 return 0; 211 } 212