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