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