1 #ifndef lint 2 static char vcid[] = "$Id: convert.c,v 1.19 1995/06/08 03:10:22 bsmith Exp bsmith $"; 3 #endif 4 5 /* Matrix conversion routines. For now, this supports only AIJ */ 6 7 #include "mpiaij.h" 8 #define ABS(a) ((a > 0) ? a : -a) 9 10 /* Determines the block diagonals within a subset of a matrix */ 11 /* For now this is just sequential -- not parallel */ 12 13 /* 14 MatDetermineDiagonals_Private - Determines the diagonal structure 15 of a matrix. 16 17 Input Parameters: 18 . mat - the matrix 19 . nb - block size 20 . irows - rows to use 21 . icols - columns to use 22 23 Output Parameters: 24 . ndiag - number of diagonals 25 . diagonals - the diagonal numbers 26 27 Note: The user must free the diagonals array. 28 */ 29 30 int MatDetermineDiagonals_Private(Mat mat,int nb,int newr,int newc, 31 int *rowrange, int *colrange,int *ndiag, int **diagonals) 32 { 33 int nd, clast, cfirst, ierr, nnc, maxd, nz, *col, *cwork, *diag; 34 int i, j, k, jdiag, cshift, row, dnew, temp; 35 Scalar *v; 36 37 VALIDHEADER(mat,MAT_COOKIE); 38 if ((newr%nb) || (newc%nb)) 39 SETERRQ(1,"MatDetermineDiagonals_Private:Invalid block size"); 40 cfirst = colrange[0]; 41 clast = colrange[newc-1]; 42 nnc = clast - cfirst + 1; 43 cwork = (int *) PETSCMALLOC( nnc * sizeof(int) ); CHKPTRQ(cwork); 44 for (i=0; i<nnc; i++) cwork[i] = -1; 45 for (i=0; i<newc; i++) cwork[colrange[i]-cfirst] = i; 46 47 /* Determine which diagonals exist: compute nd, diag[]: */ 48 /* Temporarily ssume diag[0] = 0 (main diagonal) */ 49 maxd = newr + newc - 1; /* maximum possible diagonals */ 50 diag = (int *)PETSCMALLOC( maxd * sizeof(int) ); CHKPTRQ(diag); 51 nd = 1; 52 for (i=0; i<maxd; i++) diag[i] = 0; 53 for (i=0; i<newr; i++) { 54 ierr = MatGetRow( mat, rowrange[i], &nz, &col, &v ); CHKERRQ(ierr); 55 row = i; 56 j = 0; 57 /* Skip values until we reach the first column */ 58 while (j < nz && col[j] < cfirst) j++; 59 while (j < nz) { 60 if (clast < col[j]) break; 61 cshift = cwork[col[j] - cfirst]; 62 if (cshift >= 0) { 63 /* Determine if diagonal block already exits for valid colum */ 64 dnew = 1; 65 jdiag = row/nb - cshift/nb; 66 for (k=0; k<nd; k++) { 67 if (diag[k] == jdiag) { /* diagonal exists */ 68 dnew = 0; break; 69 } 70 } 71 if (dnew) { 72 diag[nd] = jdiag; 73 nd++; 74 if (ABS(jdiag) > newr/nb) 75 { printf("ERROR jdiag\n"); } 76 } 77 } 78 j++; 79 } 80 ierr = MatRestoreRow( mat, rowrange[i], &nz, &col, &v ); CHKERRQ(ierr); 81 } 82 /* Sort diagonals in decreasing order. */ 83 for (k=0; k<nd; k++) { 84 for (j=k+1; j<nd; j++) { 85 if (diag[k] < diag[j]) { 86 temp = diag[k]; 87 diag[k] = diag[j]; 88 diag[j] = temp; 89 } 90 } 91 } 92 PETSCFREE( cwork ); 93 *ndiag = nd; 94 *diagonals = diag; 95 return 0; 96 } 97 98 /* 99 MatConvert_AIJ - Converts from MATAIJ format to another format. For 100 parallel formats, the new matrix distribution is determined by PETSc. 101 */ 102 int MatConvert_AIJ(Mat mat, MatType newtype, Mat *newmat) 103 { 104 Mat_AIJ *aij = (Mat_AIJ *) mat->data; 105 Scalar *vwork; 106 int i, ierr, nz, m = aij->m, n = aij->n, *cwork, rstart, rend; 107 108 switch (newtype) { 109 case MATROW: 110 ierr = MatCreateSequentialRow(mat->comm,m,n,0,aij->ilen,newmat); 111 CHKERRQ(ierr); break; 112 case MATMPIROW: 113 ierr = MatCreateMPIRow(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 114 m,n,0,0,0,0,newmat); /* Could do smarter memory allocation */ 115 CHKERRQ(ierr); break; 116 case MATMPIAIJ: 117 ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 118 m,n,0,0,0,0,newmat); /* Could do smarter memory allocation */ 119 CHKERRQ(ierr); break; 120 case MATDENSE: 121 ierr = MatCreateSequentialDense(mat->comm,m,n,newmat); 122 CHKERRQ(ierr); break; 123 case MATBDIAG: 124 { int nb = 1; /* Default block size = 1 */ 125 int ndiag, *diag, *rr, *cr; 126 rr = (int *) PETSCMALLOC( (m+n) * sizeof(int) ); CHKPTRQ(rr); 127 cr = rr + m; 128 for (i=0; i<m; i++) rr[i] = i; 129 for (i=0; i<n; i++) cr[i] = i; 130 OptionsGetInt(0,"-mat_bdiag_bsize",&nb); 131 ierr = MatDetermineDiagonals_Private(mat,nb,m,n,rr,cr,&ndiag,&diag); 132 CHKERRQ(ierr); 133 ierr = MatCreateSequentialBDiag(mat->comm,m,n,ndiag,nb,diag,0,newmat); 134 CHKERRQ(ierr); 135 PETSCFREE(rr), PETSCFREE(diag); 136 break; 137 } 138 case MATMPIBDIAG: 139 { int nb = 1; /* Default block size = 1 */ 140 int ndiag, *diag, *rr, *cr; 141 rr = (int *) PETSCMALLOC( (m+n) * sizeof(int) ); CHKPTRQ(rr); 142 cr = rr + m; 143 for (i=0; i<m; i++) rr[i] = i; 144 for (i=0; i<n; i++) cr[i] = i; 145 OptionsGetInt(0,"-mat_bdiag_bsize",&nb); 146 ierr = MatDetermineDiagonals_Private(mat,nb,m,n,rr,cr,&ndiag,&diag); 147 CHKERRQ(ierr); 148 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,ndiag,nb, 149 diag,0,newmat); CHKERRQ(ierr); 150 PETSCFREE(rr), PETSCFREE(diag); 151 CHKERRQ(ierr); break; 152 } 153 default: 154 SETERRQ(1,"MatConvert_AIJ:Matrix type is not currently supported"); 155 } 156 ierr = MatGetOwnershipRange(*newmat,&rstart,&rend); CHKERRQ(ierr); 157 for (i=rstart; i<rend; i++) { 158 ierr = MatGetRow(mat,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 159 ierr = MatSetValues(*newmat,1,&i,nz,cwork,vwork,INSERTVALUES); 160 CHKERRQ(ierr); 161 ierr = MatRestoreRow(mat,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 162 } 163 ierr = MatAssemblyBegin(*newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 164 ierr = MatAssemblyEnd(*newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 165 return 0; 166 } 167 /* ------------------------------------------------------------------ */ 168 /* 169 MatConvert_MPIAIJ - Converts from MATMPIAIJ format to another 170 parallel format. 171 */ 172 int MatConvert_MPIAIJ(Mat mat, MatType newtype, Mat *newmat) 173 { 174 Mat_MPIAIJ *aij = (Mat_MPIAIJ *) mat->data; 175 Mat_AIJ *Ad = (Mat_AIJ *)(aij->A->data), *Bd = (Mat_AIJ *)(aij->B->data); 176 int ierr, nz, i, ig,rstart = aij->rstart, m = aij->m, *cwork; 177 Scalar *vwork; 178 179 switch (newtype) { 180 case MATMPIROW: 181 for (i=0; i<m; i++) 182 {ierr = MatCreateMPIRow(mat->comm,m,aij->n,aij->M,aij->N,0,Ad->ilen, 183 0,Bd->ilen,newmat); CHKERRQ(ierr); } 184 break; 185 default: 186 SETERRQ(1,"MatConvert_MPIAIJ:Only MATMPIROW is currently suported"); 187 } 188 /* Each processor converts its local rows */ 189 for (i=0; i<m; i++) { 190 ig = i + rstart; 191 ierr = MatGetRow(mat,ig,&nz,&cwork,&vwork); CHKERRQ(ierr); 192 ierr = MatSetValues(*newmat,1,&ig,nz,cwork,vwork, 193 INSERTVALUES); CHKERRQ(ierr); 194 ierr = MatRestoreRow(mat,ig,&nz,&cwork,&vwork); CHKERRQ(ierr); 195 } 196 ierr = MatAssemblyBegin(*newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 197 ierr = MatAssemblyEnd(*newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 198 return 0; 199 } 200