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