1 #ifndef lint 2 static char vcid[] = "$Id: convert.c,v 1.10 1995/05/02 23:24:22 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 VALIDHEADER(mat,MAT_COOKIE); 37 if ((newr%nb) || (newc%nb)) SETERR(1,"Invalid block size."); 38 cfirst = colrange[0]; 39 clast = colrange[newc-1]; 40 nnc = clast - cfirst + 1; 41 cwork = (int *) MALLOC( nnc * sizeof(int) ); CHKPTR(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 *)MALLOC( maxd * sizeof(int) ); CHKPTR(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 ); CHKERR(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 if (abs(jdiag) > newr/nb) 73 { printf("ERROR jdiag\n"); } 74 } 75 } 76 j++; 77 } 78 ierr = MatRestoreRow( mat, rowrange[i], &nz, &col, &v ); CHKERR(ierr); 79 } 80 /* Sort diagonals in decreasing order. */ 81 for (k=0; k<nd; k++) { 82 for (j=k+1; j<nd; j++) { 83 if (diag[k] < diag[j]) { 84 temp = diag[k]; 85 diag[k] = diag[j]; 86 diag[j] = temp; 87 } 88 } 89 } 90 FREE( cwork ); 91 *ndiag = nd; 92 *diagonals = diag; 93 return 0; 94 } 95 96 /* 97 MatConvert_AIJ - Converts from MATAIJ format to another sequential format. 98 */ 99 int MatConvert_AIJ(Mat mat, MatType newtype, Mat *newmat) 100 { 101 Mat_AIJ *aij = (Mat_AIJ *) mat->data; 102 Scalar *vwork; 103 int i, ierr, nz, m = aij->m, n = aij->n, *cwork; 104 105 if (mat->type != MATAIJ) SETERR(1,"Input matrix must be MATAIJ."); 106 switch (newtype) { 107 case MATROW: 108 ierr = MatCreateSequentialRow(mat->comm,m,n,0,aij->ilen,newmat); 109 CHKERR(ierr); break; 110 case MATDENSE: 111 ierr = MatCreateSequentialDense(mat->comm,m,n,newmat); 112 CHKERR(ierr); break; 113 case MATBDIAG: 114 { int nb = 1; /* Default block size = 1 */ 115 int ndiag, *diag, *rr, *cr; 116 rr = (int *) MALLOC( (m+n) * sizeof(int) ); CHKPTR(rr); 117 cr = rr + m; 118 for (i=0; i<m; i++) rr[i] = i; 119 for (i=0; i<n; i++) cr[i] = i; 120 OptionsGetInt(0,0,"-mat_bdiag_bsize",&nb); 121 ierr = MatDetermineDiagonals_Private(mat,nb,m,n,rr,cr,&ndiag,&diag); 122 ierr = MatCreateSequentialBDiag(mat->comm,m,n,ndiag,nb,diag,0,newmat); 123 124 MatAssemblyEnd(*newmat,FINAL_ASSEMBLY); MatView(*newmat,0); 125 126 FREE(rr), FREE(diag); 127 CHKERR(ierr); break; 128 } 129 default: 130 SETERR(1,"Only MATROW, MATDENSE, and MATBDIAG are currently supported."); 131 } 132 /* for (i=0; i<m; i++) { */ 133 for (i=0; i<1; i++) { 134 ierr = MatGetRow(mat,i,&nz,&cwork,&vwork); CHKERR(ierr); 135 ierr = MatSetValues(*newmat,1,&i,nz,cwork,vwork,INSERTVALUES); 136 CHKERR(ierr); 137 ierr = MatRestoreRow(mat,i,&nz,&cwork,&vwork); CHKERR(ierr); 138 } 139 ierr = MatAssemblyBegin(*newmat,FINAL_ASSEMBLY); CHKERR(ierr); 140 ierr = MatAssemblyEnd(*newmat,FINAL_ASSEMBLY); CHKERR(ierr); 141 return 0; 142 } 143 /* ------------------------------------------------------------------ */ 144 /* 145 MatConvert_MPIAIJ - Converts from MATMPIAIJ format to another 146 parallel format. 147 */ 148 int MatConvert_MPIAIJ(Mat mat, MatType newtype, Mat *newmat) 149 { 150 Mat_MPIAIJ *aij = (Mat_MPIAIJ *) mat->data; 151 Mat_AIJ *Ad = (Mat_AIJ *)(aij->A->data), *Bd = (Mat_AIJ *)(aij->B->data); 152 int ierr, nz, i, ig,rstart = aij->rstart, m = aij->m, *cwork; 153 Scalar *vwork; 154 155 if (mat->type != MATMPIAIJ) SETERR(1,"Input matrix must be MATMPIAIJ."); 156 switch (newtype) { 157 case MATMPIROW: 158 for (i=0; i<m; i++) 159 {ierr = MatCreateMPIRow(mat->comm,m,aij->n,aij->M,aij->N,0,Ad->ilen, 160 0,Bd->ilen,newmat); CHKERR(ierr); } 161 break; 162 default: 163 SETERR(1,"Only MATMPIROW is currently suported."); 164 } 165 /* Each processor converts its local rows */ 166 for (i=0; i<m; i++) { 167 ig = i + rstart; 168 ierr = MatGetRow(mat,ig,&nz,&cwork,&vwork); CHKERR(ierr); 169 ierr = MatSetValues(*newmat,1,&ig,nz,cwork,vwork, 170 INSERTVALUES); CHKERR(ierr); 171 ierr = MatRestoreRow(mat,ig,&nz,&cwork,&vwork); CHKERR(ierr); 172 } 173 ierr = MatAssemblyBegin(*newmat,FINAL_ASSEMBLY); CHKERR(ierr); 174 ierr = MatAssemblyEnd(*newmat,FINAL_ASSEMBLY); CHKERR(ierr); 175 return 0; 176 } 177