1 2 #include "src/mat/impls/aij/seq/aij.h" 3 #include "src/mat/impls/baij/seq/baij.h" 4 #include "src/mat/impls/sbaij/seq/sbaij.h" 5 6 EXTERN_C_BEGIN 7 #undef __FUNCT__ 8 #define __FUNCT__ "MatConvert_SeqSBAI_SeqAIJ" 9 PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A,const MatType newtype,MatReuse reuse,Mat *newmat) 10 { 11 Mat B; 12 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 13 Mat_SeqAIJ *b; 14 PetscErrorCode ierr; 15 PetscInt *ai=a->i,*aj=a->j,m=A->m,n=A->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp; 16 PetscInt bs=A->bs,bs2=bs*bs,mbs=A->m/bs; 17 PetscScalar *av,*bv; 18 19 PetscFunctionBegin; 20 21 /* compute rowlengths of newmat */ 22 ierr = PetscMalloc((2*m+1)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 23 rowstart = rowlengths + m; 24 25 for (i=0; i<mbs; i++) rowlengths[i*bs] = 0; 26 aj = a->j; 27 k = 0; 28 for (i=0; i<mbs; i++) { 29 nz = ai[i+1] - ai[i]; 30 aj++; /* skip diagonal */ 31 for (j=1; j<nz; j++) { /* no. of lower triangular blocks */ 32 rowlengths[(*aj)*bs]++; aj++; 33 } 34 rowlengths[k] += nz; /* no. of upper triangular blocks */ 35 rowlengths[k] *= bs; 36 for (j=1; j<bs; j++) { 37 rowlengths[k+j] = rowlengths[k]; 38 } 39 k += bs; 40 /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */ 41 } 42 43 ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr); 44 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 45 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 46 ierr = MatSetOption(B,MAT_COLUMN_ORIENTED);CHKERRQ(ierr); 47 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 48 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 49 B->bs = A->bs; 50 51 b = (Mat_SeqAIJ*)(B->data); 52 bi = b->i; 53 bj = b->j; 54 bv = b->a; 55 56 /* set b->i */ 57 bi[0] = 0; rowstart[0] = 0; 58 for (i=0; i<mbs; i++){ 59 for (j=0; j<bs; j++){ 60 b->ilen[i*bs+j] = rowlengths[i*bs]; 61 rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs]; 62 } 63 bi[i+1] = bi[i] + rowlengths[i*bs]/bs; 64 } 65 if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz-mbs: %D\n",bi[mbs],2*a->nz - mbs); 66 67 /* set b->j and b->a */ 68 aj = a->j; av = a->a; 69 for (i=0; i<mbs; i++) { 70 /* diagonal block */ 71 for (j=0; j<bs; j++){ /* row i*bs+j */ 72 itmp = i*bs+j; 73 for (k=0; k<bs; k++){ /* col i*bs+k */ 74 *(bj + rowstart[itmp]) = (*aj)*bs+k; 75 *(bv + rowstart[itmp]) = *(av+k*bs+j); 76 rowstart[itmp]++; 77 } 78 } 79 aj++; av += bs2; 80 81 nz = ai[i+1] - ai[i] -1; 82 while (nz--){ 83 /* lower triangular blocks */ 84 for (j=0; j<bs; j++){ /* row (*aj)*bs+j */ 85 itmp = (*aj)*bs+j; 86 for (k=0; k<bs; k++){ /* col i*bs+k */ 87 *(bj + rowstart[itmp]) = i*bs+k; 88 *(bv + rowstart[itmp]) = *(av+k*bs+j); 89 rowstart[itmp]++; 90 } 91 } 92 /* upper triangular blocks */ 93 for (j=0; j<bs; j++){ /* row i*bs+j */ 94 itmp = i*bs+j; 95 for (k=0; k<bs; k++){ /* col (*aj)*bs+k */ 96 *(bj + rowstart[itmp]) = (*aj)*bs+k; 97 *(bv + rowstart[itmp]) = *(av+k*bs+j); 98 rowstart[itmp]++; 99 } 100 } 101 aj++; av += bs2; 102 } 103 } 104 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 105 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 106 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 107 108 if (reuse == MAT_REUSE_MATRIX) { 109 ierr = MatHeaderCopy(A,B);CHKERRQ(ierr); 110 } else { 111 *newmat = B; 112 } 113 PetscFunctionReturn(0); 114 } 115 #undef __FUNCT__ 116 #define __FUNCT__ "MatConvert_SeqAIJ_SeqSBAIJ" 117 PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,const MatType newtype,MatReuse reuse,Mat *newmat) { 118 Mat B; 119 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 120 Mat_SeqSBAIJ *b; 121 PetscErrorCode ierr; 122 PetscInt *ai=a->i,*aj,m=A->M,n=A->N,i,j,*bi,*bj,*rowlengths; 123 PetscScalar *av,*bv; 124 125 PetscFunctionBegin; 126 if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 127 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */ 128 129 ierr = PetscMalloc(m*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 130 for (i=0; i<m; i++) { 131 rowlengths[i] = ai[i+1] - a->diag[i]; 132 } 133 ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr); 134 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 135 ierr = MatSeqSBAIJSetPreallocation(B,1,0,rowlengths);CHKERRQ(ierr); 136 137 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 138 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 139 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 140 141 b = (Mat_SeqSBAIJ*)(B->data); 142 bi = b->i; 143 bj = b->j; 144 bv = b->a; 145 146 bi[0] = 0; 147 for (i=0; i<m; i++) { 148 aj = a->j + a->diag[i]; 149 av = a->a + a->diag[i]; 150 for (j=0; j<rowlengths[i]; j++){ 151 *bj = *aj; bj++; aj++; 152 *bv = *av; bv++; av++; 153 } 154 bi[i+1] = bi[i] + rowlengths[i]; 155 b->ilen[i] = rowlengths[i]; 156 } 157 158 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 159 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 160 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 161 162 if (reuse == MAT_REUSE_MATRIX) { 163 ierr = MatHeaderCopy(A,B);CHKERRQ(ierr); 164 } else { 165 *newmat = B; 166 } 167 PetscFunctionReturn(0); 168 } 169 EXTERN_C_END 170 171 EXTERN_C_BEGIN 172 #undef __FUNCT__ 173 #define __FUNCT__ "MatConvert_SeqSBAI_SeqBAIJ" 174 PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A,const MatType newtype,MatReuse reuse,Mat *newmat) 175 { 176 Mat B; 177 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 178 Mat_SeqBAIJ *b; 179 PetscErrorCode ierr; 180 PetscInt *ai=a->i,*aj=a->j,m=A->m,n=A->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp; 181 PetscInt bs=A->bs,bs2=bs*bs,mbs=m/bs; 182 PetscScalar *av,*bv; 183 184 PetscFunctionBegin; 185 /* compute browlengths of newmat */ 186 ierr = PetscMalloc(2*mbs*sizeof(PetscInt),&browlengths);CHKERRQ(ierr); 187 browstart = browlengths + mbs; 188 for (i=0; i<mbs; i++) browlengths[i] = 0; 189 aj = a->j; 190 for (i=0; i<mbs; i++) { 191 nz = ai[i+1] - ai[i]; 192 aj++; /* skip diagonal */ 193 for (k=1; k<nz; k++) { /* no. of lower triangular blocks */ 194 browlengths[*aj]++; aj++; 195 } 196 browlengths[i] += nz; /* no. of upper triangular blocks */ 197 } 198 199 ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr); 200 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 201 ierr = MatSeqBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 202 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 203 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 204 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 205 206 b = (Mat_SeqBAIJ*)(B->data); 207 bi = b->i; 208 bj = b->j; 209 bv = b->a; 210 211 /* set b->i */ 212 bi[0] = 0; 213 for (i=0; i<mbs; i++){ 214 b->ilen[i] = browlengths[i]; 215 bi[i+1] = bi[i] + browlengths[i]; 216 browstart[i] = bi[i]; 217 } 218 if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz - mbs: %D\n",bi[mbs],2*a->nz - mbs); 219 220 /* set b->j and b->a */ 221 aj = a->j; av = a->a; 222 for (i=0; i<mbs; i++) { 223 /* diagonal block */ 224 *(bj + browstart[i]) = *aj; aj++; 225 itmp = bs2*browstart[i]; 226 for (k=0; k<bs2; k++){ 227 *(bv + itmp + k) = *av; av++; 228 } 229 browstart[i]++; 230 231 nz = ai[i+1] - ai[i] -1; 232 while (nz--){ 233 /* lower triangular blocks */ 234 *(bj + browstart[*aj]) = i; 235 itmp = bs2*browstart[*aj]; 236 for (k=0; k<bs2; k++){ 237 *(bv + itmp + k) = *(av + k); 238 } 239 browstart[*aj]++; 240 241 /* upper triangular blocks */ 242 *(bj + browstart[i]) = *aj; aj++; 243 itmp = bs2*browstart[i]; 244 for (k=0; k<bs2; k++){ 245 *(bv + itmp + k) = *av; av++; 246 } 247 browstart[i]++; 248 } 249 } 250 ierr = PetscFree(browlengths);CHKERRQ(ierr); 251 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 252 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 253 254 if (reuse == MAT_REUSE_MATRIX) { 255 ierr = MatHeaderCopy(A,B);CHKERRQ(ierr); 256 } else { 257 *newmat = B; 258 } 259 PetscFunctionReturn(0); 260 } 261 #undef __FUNCT__ 262 #define __FUNCT__ "MatConvert_SeqBAIJ_SeqSBAIJ" 263 PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A,const MatType newtype,MatReuse reuse,Mat *newmat) 264 { 265 Mat B; 266 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 267 Mat_SeqSBAIJ *b; 268 PetscErrorCode ierr; 269 PetscInt *ai=a->i,*aj,m=A->m,n=A->n,i,j,k,*bi,*bj,*browlengths; 270 PetscInt bs=A->bs,bs2=bs*bs,mbs=m/bs; 271 PetscScalar *av,*bv; 272 273 PetscFunctionBegin; 274 if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 275 ierr = MatMissingDiagonal_SeqBAIJ(A);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */ 276 277 ierr = PetscMalloc(mbs*sizeof(PetscInt),&browlengths);CHKERRQ(ierr); 278 for (i=0; i<mbs; i++) { 279 browlengths[i] = ai[i+1] - a->diag[i]; 280 } 281 282 ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr); 283 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 284 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 285 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 286 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 287 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 288 289 b = (Mat_SeqSBAIJ*)(B->data); 290 bi = b->i; 291 bj = b->j; 292 bv = b->a; 293 294 bi[0] = 0; 295 for (i=0; i<mbs; i++) { 296 aj = a->j + a->diag[i]; 297 av = a->a + (a->diag[i])*bs2; 298 for (j=0; j<browlengths[i]; j++){ 299 *bj = *aj; bj++; aj++; 300 for (k=0; k<bs2; k++){ 301 *bv = *av; bv++; av++; 302 } 303 } 304 bi[i+1] = bi[i] + browlengths[i]; 305 b->ilen[i] = browlengths[i]; 306 } 307 ierr = PetscFree(browlengths);CHKERRQ(ierr); 308 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 309 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 310 311 if (reuse == MAT_REUSE_MATRIX) { 312 ierr = MatHeaderCopy(A);CHKERRQ(ierr); 313 } else { 314 *newmat = B; 315 } 316 317 PetscFunctionReturn(0); 318 } 319 EXTERN_C_END 320