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