1 #ifndef lint 2 static char vcid[] = "$Id: convert.c,v 1.41 1996/01/12 22:08:29 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 MATSEQROW: 27 ierr = MatCreateSeqRow(mat->comm,m,n,0,PETSC_NULL,M); CHKERRQ(ierr); 28 break; 29 case MATMPIROW: 30 ierr = MatCreateMPIRow(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 31 m,n,0,PETSC_NULL,0,PETSC_NULL,M); CHKERRQ(ierr); 32 break; 33 case MATMPIROWBS: 34 if (m != n) SETERRQ(1,"MatConvert:MATMPIROWBS matrix must be square"); 35 ierr = MatCreateMPIRowbs(MPI_COMM_WORLD,PETSC_DECIDE,m,0,PETSC_NULL, 36 PETSC_NULL,M); CHKERRQ(ierr); 37 break; 38 case MATMPIAIJ: 39 ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 40 m,n,0,PETSC_NULL,0,PETSC_NULL,M); CHKERRQ(ierr); 41 break; 42 case MATSEQDENSE: 43 ierr = MatCreateSeqDense(mat->comm,m,n,PETSC_NULL,M); CHKERRQ(ierr); 44 break; 45 case MATMPIDENSE: 46 ierr = MatCreateMPIDense(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 47 m,n,PETSC_NULL,M); CHKERRQ(ierr); 48 break; 49 case MATSEQBDIAG: 50 { 51 int nb = 1; /* Default block size = 1 */ 52 ierr = OptionsGetInt(PETSC_NULL,"-mat_bdiag_bsize",&nb,&flg); CHKERRQ(ierr); 53 ierr = MatCreateSeqBDiag(mat->comm,m,n,0,nb,PETSC_NULL,PETSC_NULL,M); CHKERRQ(ierr); 54 break; 55 } 56 case MATMPIBDIAG: 57 { 58 int nb = 1; /* Default block size = 1 */ 59 ierr = OptionsGetInt(PETSC_NULL,"-mat_bdiag_bsize",&nb,&flg); CHKERRQ(ierr); 60 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,0,nb,PETSC_NULL, 61 PETSC_NULL,M); CHKERRQ(ierr); 62 break; 63 } 64 default: 65 SETERRQ(1,"MatConvert:Matrix type is not currently supported"); 66 } 67 ierr = MatGetOwnershipRange(*M,&rstart,&rend); CHKERRQ(ierr); 68 for (i=rstart; i<rend; i++) { 69 ierr = MatGetRow(mat,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 70 ierr = MatSetValues(*M,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 71 ierr = MatRestoreRow(mat,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 72 } 73 ierr = MatAssemblyBegin(*M,FINAL_ASSEMBLY); CHKERRQ(ierr); 74 ierr = MatAssemblyEnd(*M,FINAL_ASSEMBLY); CHKERRQ(ierr); 75 return 0; 76 } 77 /* -------------------------------------------------------------- */ 78 /* 79 MatConvert_SeqAIJ - Converts from MATSEQAIJ format to another format. For 80 parallel formats, the new matrix distribution is determined by PETSc. 81 */ 82 int MatConvert_SeqAIJ(Mat A, MatType newtype, Mat *B) 83 { 84 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 85 Scalar *vwork; 86 int i, ierr, nz, m = a->m, n = a->n, *cwork, rstart, rend,flg; 87 88 switch (newtype) { 89 case MATSEQROW: 90 ierr = MatCreateSeqRow(A->comm,m,n,0,a->ilen,B); CHKERRQ(ierr); 91 break; 92 case MATMPIROW: 93 ierr = MatCreateMPIRow(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 94 m,n,0,PETSC_NULL,0,PETSC_NULL,B); CHKERRQ(ierr); 95 break; 96 case MATMPIROWBS: 97 if (m != n) SETERRQ(1,"MatConvert_SeqAIJ:MATMPIROWBS matrix must be square"); 98 ierr = MatCreateMPIRowbs(MPI_COMM_WORLD,PETSC_DECIDE,m,0,PETSC_NULL,PETSC_NULL,B); 99 CHKERRQ(ierr); 100 break; 101 case MATMPIAIJ: 102 ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 103 m,n,0,PETSC_NULL,0,PETSC_NULL,B); CHKERRQ(ierr); 104 break; 105 case MATSEQDENSE: 106 ierr = MatCreateSeqDense(A->comm,m,n,PETSC_NULL,B); CHKERRQ(ierr); 107 break; 108 case MATMPIDENSE: 109 ierr = MatCreateMPIDense(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 110 m,n,PETSC_NULL,B); CHKERRQ(ierr); 111 break; 112 case MATSEQBDIAG: 113 { 114 int nb = 1; /* Default block size = 1 */ 115 ierr = OptionsGetInt(PETSC_NULL,"-mat_bdiag_bsize",&nb,&flg); CHKERRQ(ierr); 116 ierr = MatCreateSeqBDiag(A->comm,m,n,0,nb,PETSC_NULL,PETSC_NULL,B); CHKERRQ(ierr); 117 break; 118 } 119 case MATMPIBDIAG: 120 { 121 int nb = 1; /* Default block size = 1 */ 122 ierr = OptionsGetInt(PETSC_NULL,"-mat_bdiag_bsize",&nb,&flg); CHKERRQ(ierr); 123 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,0,nb,PETSC_NULL, 124 PETSC_NULL,B); CHKERRQ(ierr); 125 break; 126 } 127 default: 128 SETERRQ(1,"MatConvert_SeqAIJ:Matrix type is not currently supported"); 129 } 130 ierr = MatGetOwnershipRange(*B,&rstart,&rend); CHKERRQ(ierr); 131 for (i=rstart; i<rend; i++) { 132 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 133 ierr = MatSetValues(*B,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 134 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 135 } 136 ierr = MatAssemblyBegin(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 137 ierr = MatAssemblyEnd(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 138 return 0; 139 } 140 /* ------------------------------------------------------------------ */ 141 /* 142 MatConvert_MPIAIJ - Converts from MATMPIAIJ format to another 143 parallel format. 144 */ 145 int MatConvert_MPIAIJ(Mat A, MatType newtype, Mat *B) 146 { 147 Mat_MPIAIJ *a = (Mat_MPIAIJ *) A->data; 148 Mat_SeqAIJ *Ad = (Mat_SeqAIJ *)(a->A->data), *Bd = (Mat_SeqAIJ *)(a->B->data); 149 int ierr, nz, i, ig, rstart = a->rstart, m = a->m, *cwork; 150 Scalar *vwork; 151 152 switch (newtype) { 153 case MATMPIROW: 154 ierr = MatCreateMPIRow(A->comm,m,a->n,a->M,a->N,0,Ad->ilen, 155 0,Bd->ilen,B); CHKERRQ(ierr); 156 break; 157 default: 158 SETERRQ(1,"MatConvert_MPIAIJ:Only MATMPIROW is currently suported"); 159 } 160 /* Each processor converts its local rows */ 161 for (i=0; i<m; i++) { 162 ig = i + rstart; 163 ierr = MatGetRow(A,ig,&nz,&cwork,&vwork); CHKERRQ(ierr); 164 ierr = MatSetValues(*B,1,&ig,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 165 ierr = MatRestoreRow(A,ig,&nz,&cwork,&vwork); CHKERRQ(ierr); 166 } 167 ierr = MatAssemblyBegin(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 168 ierr = MatAssemblyEnd(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 169 return 0; 170 } 171 /* ------------------------------------------------------------------ */ 172 /* 173 MatConvert_SeqBDiag - Converts from MATSEQBDiag format to another format. For 174 parallel formats, the new matrix distribution is determined by PETSc. 175 */ 176 int MatConvert_SeqBDiag(Mat A, MatType newtype, Mat *B) 177 { 178 Mat_SeqBDiag *a = (Mat_SeqBDiag *) A->data; 179 Scalar *vwork, *vw2; 180 int i, ierr, nz, m = a->m, n = a->n, *cwork, rstart, rend; 181 int j, *cw2, ict; 182 183 /* rough over-estimate; could refine for individual rows */ 184 nz = PetscMin(n,a->nd*a->nb); 185 switch (newtype) { 186 case MATSEQAIJ: 187 ierr = MatCreateSeqAIJ(A->comm,m,n,nz,PETSC_NULL,B); CHKERRQ(ierr); 188 break; 189 case MATSEQROW: 190 ierr = MatCreateSeqRow(A->comm,m,n,nz,PETSC_NULL,B); CHKERRQ(ierr); 191 break; 192 case MATMPIROW: 193 ierr = MatCreateMPIRow(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 194 m,n,0,PETSC_NULL,0,PETSC_NULL,B); CHKERRQ(ierr); 195 break; 196 case MATMPIROWBS: 197 if (m != n) SETERRQ(1,"MatConvert_SeqBDiag:MATMPIROWBS matrix must be square"); 198 ierr = MatCreateMPIRowbs(MPI_COMM_WORLD,PETSC_DECIDE,m,0,PETSC_NULL,PETSC_NULL, 199 B); CHKERRQ(ierr); 200 break; 201 case MATMPIAIJ: 202 ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 203 m,n,0,PETSC_NULL,0,PETSC_NULL,B); CHKERRQ(ierr); 204 break; 205 case MATSEQDENSE: 206 ierr = MatCreateSeqDense(A->comm,m,n,PETSC_NULL,B); CHKERRQ(ierr); 207 break; 208 case MATMPIDENSE: 209 ierr = MatCreateMPIDense(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 210 m,n,PETSC_NULL,B); CHKERRQ(ierr); 211 break; 212 case MATMPIBDIAG: 213 { 214 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,a->nd,a->nb, 215 PETSC_NULL,PETSC_NULL,B); CHKERRQ(ierr); 216 break; 217 } 218 default: 219 SETERRQ(1,"MatConvert_SeqBDiag:Matrix type is not currently supported"); 220 } 221 ierr = MatGetOwnershipRange(*B,&rstart,&rend); CHKERRQ(ierr); 222 223 cw2 = (int *)PetscMalloc( n * sizeof(int) ); CHKPTRQ(cw2); 224 vw2 = (Scalar *)PetscMalloc( n * sizeof(Scalar) ); CHKPTRQ(vw2); 225 for (i=rstart; i<rend; i++) { 226 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 227 ict = 0; /* strip out the zero elements ... is this what we really want? */ 228 for (j=0; j<nz; j++) { 229 if (vwork[j] != 0) {vw2[ict] = vwork[j]; cw2[ict] = cwork[j]; ict++;} 230 } 231 if (ict) 232 {ierr = MatSetValues(*B,1,&i,ict,cw2,vw2,INSERT_VALUES); CHKERRQ(ierr);} 233 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 234 } 235 PetscFree(cw2); PetscFree(vw2); 236 ierr = MatAssemblyBegin(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 237 ierr = MatAssemblyEnd(*B,FINAL_ASSEMBLY); CHKERRQ(ierr); 238 return 0; 239 } 240