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