1 #ifndef lint 2 static char vcid[] = "$Id: convert.c,v 1.46 1996/02/13 23:30:02 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_block_size",&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_block_size",&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_block_size",&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_block_size",&nb,&flg);CHKERRQ(ierr); 109 CHKERRQ(ierr); 110 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,0,nb,PETSC_NULL, 111 PETSC_NULL,B); CHKERRQ(ierr); 112 break; 113 } 114 default: 115 SETERRQ(1,"MatConvert_SeqAIJ:Matrix type is not currently supported"); 116 } 117 ierr = MatGetOwnershipRange(*B,&rstart,&rend); CHKERRQ(ierr); 118 for (i=rstart; i<rend; i++) { 119 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 120 ierr = MatSetValues(*B,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 121 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 122 } 123 ierr = MatAssemblyBegin(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 124 ierr = MatAssemblyEnd(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 125 return 0; 126 } 127 /* ------------------------------------------------------------------ */ 128 /* 129 MatConvert_MPIAIJ - Converts from MATMPIAIJ format to another 130 parallel format. 131 */ 132 int MatConvert_MPIAIJ(Mat A, MatType newtype, Mat *B) 133 { 134 SETERRQ(1,"MatConvert_MPIAIJ:Not currently suported"); 135 136 /* Each processor converts its local rows */ 137 /* ---------------------------------------------------- 138 Mat_MPIAIJ *a = (Mat_MPIAIJ *) A->data; 139 int ierr, nz, i, ig, rstart = a->rstart, m = a->m, *cwork; 140 Scalar *vwork; 141 142 for (i=0; i<m; i++) { 143 ig = i + rstart; 144 ierr = MatGetRow(A,ig,&nz,&cwork,&vwork); CHKERRQ(ierr); 145 ierr = MatSetValues(*B,1,&ig,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 146 ierr = MatRestoreRow(A,ig,&nz,&cwork,&vwork); CHKERRQ(ierr); 147 } 148 ierr = MatAssemblyBegin(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 149 ierr = MatAssemblyEnd(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 150 return 0; 151 ---------------------------------------------------------*/ 152 } 153 /* ------------------------------------------------------------------ */ 154 /* 155 MatConvert_SeqBDiag - Converts from MATSEQBDiag format to another format. For 156 parallel formats, the new matrix distribution is determined by PETSc. 157 */ 158 int MatConvert_SeqBDiag(Mat A, MatType newtype, Mat *B) 159 { 160 Mat_SeqBDiag *a = (Mat_SeqBDiag *) A->data; 161 Scalar *vwork, *vw2; 162 int i, ierr, nz, m = a->m, n = a->n, *cwork, rstart, rend; 163 int j, *cw2, ict; 164 165 /* rough over-estimate; could refine for individual rows */ 166 nz = PetscMin(n,a->nd*a->nb); 167 switch (newtype) { 168 case MATSEQAIJ: 169 ierr = MatCreateSeqAIJ(A->comm,m,n,nz,PETSC_NULL,B); CHKERRQ(ierr); 170 break; 171 case MATMPIROWBS: 172 if (m != n) SETERRQ(1,"MatConvert_SeqBDiag:MATMPIROWBS matrix must be square"); 173 ierr = MatCreateMPIRowbs(MPI_COMM_WORLD,PETSC_DECIDE,m,0,PETSC_NULL,PETSC_NULL, 174 B); CHKERRQ(ierr); 175 break; 176 case MATMPIAIJ: 177 ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 178 m,n,0,PETSC_NULL,0,PETSC_NULL,B); CHKERRQ(ierr); 179 break; 180 case MATSEQDENSE: 181 ierr = MatCreateSeqDense(A->comm,m,n,PETSC_NULL,B); CHKERRQ(ierr); 182 break; 183 case MATMPIDENSE: 184 ierr = MatCreateMPIDense(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 185 m,n,PETSC_NULL,B); CHKERRQ(ierr); 186 break; 187 case MATMPIBDIAG: 188 { 189 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,a->nd,a->nb, 190 PETSC_NULL,PETSC_NULL,B); CHKERRQ(ierr); 191 break; 192 } 193 default: 194 SETERRQ(1,"MatConvert_SeqBDiag:Matrix type is not currently supported"); 195 } 196 ierr = MatGetOwnershipRange(*B,&rstart,&rend); CHKERRQ(ierr); 197 198 cw2 = (int *)PetscMalloc( n * sizeof(int) ); CHKPTRQ(cw2); 199 vw2 = (Scalar *)PetscMalloc( n * sizeof(Scalar) ); CHKPTRQ(vw2); 200 for (i=rstart; i<rend; i++) { 201 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 202 ict = 0; /* strip out the zero elements ... is this what we really want? */ 203 for (j=0; j<nz; j++) { 204 if (vwork[j] != 0) {vw2[ict] = vwork[j]; cw2[ict] = cwork[j]; ict++;} 205 } 206 if (ict) { 207 ierr = MatSetValues(*B,1,&i,ict,cw2,vw2,INSERT_VALUES); CHKERRQ(ierr); 208 } 209 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 210 } 211 PetscFree(cw2); PetscFree(vw2); 212 ierr = MatAssemblyBegin(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 213 ierr = MatAssemblyEnd(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 214 return 0; 215 } 216