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