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