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