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