1 #ifndef lint 2 static char vcid[] = "$Id: convert.c,v 1.26 1995/09/21 20:11:31 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,int *rowrange, 30 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)) SETERRQ(1,"MatDetermineDiagonals_Private:Bad block size"); 38 cfirst = colrange[0]; 39 clast = colrange[newc-1]; 40 nnc = clast - cfirst + 1; 41 cwork = (int *) PETSCMALLOC( nnc * sizeof(int) ); CHKPTRQ(cwork); 42 for (i=0; i<nnc; i++) cwork[i] = -1; 43 for (i=0; i<newc; i++) cwork[colrange[i]-cfirst] = i; 44 45 /* Determine which diagonals exist: compute nd, diag[]: */ 46 /* Temporarily ssume diag[0] = 0 (main diagonal) */ 47 maxd = newr + newc - 1; /* maximum possible diagonals */ 48 diag = (int *)PETSCMALLOC( maxd * sizeof(int) ); CHKPTRQ(diag); 49 nd = 1; 50 for (i=0; i<maxd; i++) diag[i] = 0; 51 for (i=0; i<newr; i++) { 52 ierr = MatGetRow( mat, rowrange[i], &nz, &col, &v ); CHKERRQ(ierr); 53 row = i; 54 j = 0; 55 /* Skip values until we reach the first column */ 56 while (j < nz && col[j] < cfirst) j++; 57 while (j < nz) { 58 if (clast < col[j]) break; 59 cshift = cwork[col[j] - cfirst]; 60 if (cshift >= 0) { 61 /* Determine if diagonal block already exits for valid colum */ 62 dnew = 1; 63 jdiag = row/nb - cshift/nb; 64 for (k=0; k<nd; k++) { 65 if (diag[k] == jdiag) { /* diagonal exists */ 66 dnew = 0; break; 67 } 68 } 69 if (dnew) { 70 diag[nd] = jdiag; 71 nd++; 72 /* what the hell is this? Not good PETSc style */ 73 if (PETSCABS(jdiag) > newr/nb) { 74 SETERRQ(1," MatDetermineDiagonals_Private: bad jdiag"); 75 } 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_SeqAIJ - Converts from MATSEQAIJ format to another format. For 100 parallel formats, the new matrix distribution is determined by PETSc. 101 */ 102 int MatConvert_SeqAIJ(Mat mat, MatType newtype, Mat *newmat) 103 { 104 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) 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 MATSEQROW: 110 ierr = MatCreateSeqRow(mat->comm,m,n,0,aij->ilen,newmat);CHKERRQ(ierr); 111 break; 112 case MATMPIROW: 113 if (m != n) SETERRQ(1,"MatConvert_SeqAIJ: MPIRowbs matrix must be square"); 114 ierr = MatCreateMPIRow(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 115 m,n,0,0,0,0,newmat); CHKERRQ(ierr); 116 break; 117 case MATMPIROWBS: 118 ierr = MatCreateMPIRowbs(MPI_COMM_WORLD,PETSC_DECIDE, 119 m,0,0,0,newmat);CHKERRQ(ierr); 120 break; 121 case MATMPIAIJ: 122 ierr = MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, 123 m,n,0,0,0,0,newmat);CHKERRQ(ierr); 124 break; 125 case MATSEQDENSE: 126 ierr = MatCreateSeqDense(mat->comm,m,n,newmat); CHKERRQ(ierr); 127 break; 128 case MATSEQBDIAG: 129 { 130 int nb = 1, /* Default block size = 1 */ ndiag, *diag, *rr, *cr; 131 rr = (int *) PETSCMALLOC( (m+n) * sizeof(int) ); CHKPTRQ(rr); 132 cr = rr + m; 133 for (i=0; i<m; i++) rr[i] = i; 134 for (i=0; i<n; i++) cr[i] = i; 135 OptionsGetInt(0,"-mat_bdiag_bsize",&nb); 136 ierr = MatDetermineDiagonals_Private(mat,nb,m,n,rr,cr,&ndiag,&diag);CHKERRQ(ierr); 137 ierr = MatCreateSeqBDiag(mat->comm,m,n,ndiag,nb,diag,0,newmat); CHKERRQ(ierr); 138 PETSCFREE(rr), PETSCFREE(diag); 139 break; 140 } 141 case MATMPIBDIAG: 142 { 143 int nb = 1, /* Default block size = 1 */ ndiag, *diag, *rr, *cr; 144 rr = (int *) PETSCMALLOC( (m+n) * sizeof(int) ); CHKPTRQ(rr); 145 cr = rr + m; 146 for (i=0; i<m; i++) rr[i] = i; 147 for (i=0; i<n; i++) cr[i] = i; 148 OptionsGetInt(0,"-mat_bdiag_bsize",&nb); 149 ierr = MatDetermineDiagonals_Private(mat,nb,m,n,rr,cr,&ndiag,&diag); 150 CHKERRQ(ierr); 151 ierr = MatCreateMPIBDiag(MPI_COMM_WORLD,PETSC_DECIDE,m,n,ndiag,nb, 152 diag,0,newmat); CHKERRQ(ierr); 153 PETSCFREE(rr), PETSCFREE(diag); 154 break; 155 } 156 default: 157 SETERRQ(1,"MatConvert_SeqAIJ:Matrix type is not currently supported"); 158 } 159 ierr = MatGetOwnershipRange(*newmat,&rstart,&rend); CHKERRQ(ierr); 160 for (i=rstart; i<rend; i++) { 161 ierr = MatGetRow(mat,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 162 ierr = MatSetValues(*newmat,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 163 ierr = MatRestoreRow(mat,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 164 } 165 ierr = MatAssemblyBegin(*newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 166 ierr = MatAssemblyEnd(*newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 167 return 0; 168 } 169 /* ------------------------------------------------------------------ */ 170 /* 171 MatConvert_MPIAIJ - Converts from MATMPIAIJ format to another 172 parallel format. 173 */ 174 int MatConvert_MPIAIJ(Mat mat, MatType newtype, Mat *newmat) 175 { 176 Mat_MPIAIJ *aij = (Mat_MPIAIJ *) mat->data; 177 Mat_SeqAIJ *Ad = (Mat_SeqAIJ *)(aij->A->data), *Bd = (Mat_SeqAIJ *)(aij->B->data); 178 int ierr, nz, i, ig,rstart = aij->rstart, m = aij->m, *cwork; 179 Scalar *vwork; 180 181 switch (newtype) { 182 case MATMPIROW: 183 ierr = MatCreateMPIRow(mat->comm,m,aij->n,aij->M,aij->N,0,Ad->ilen, 184 0,Bd->ilen,newmat); CHKERRQ(ierr); 185 break; 186 default: 187 SETERRQ(1,"MatConvert_MPIAIJ:Only MATMPIROW is currently suported"); 188 } 189 /* Each processor converts its local rows */ 190 for (i=0; i<m; i++) { 191 ig = i + rstart; 192 ierr = MatGetRow(mat,ig,&nz,&cwork,&vwork); CHKERRQ(ierr); 193 ierr = MatSetValues(*newmat,1,&ig,nz,cwork,vwork,INSERT_VALUES); 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 201 202