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_SeqSBAIJ_SeqAIJ" 10 PetscErrorCode 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,diagcnt=0; 18 MatScalar *av,*bv; 19 20 PetscFunctionBegin; 21 /* compute rowlengths of newmat */ 22 ierr = PetscMalloc2(m,PetscInt,&rowlengths,m+1,PetscInt,&rowstart);CHKERRQ(ierr); 23 24 for (i=0; i<mbs; i++) rowlengths[i*bs] = 0; 25 aj = a->j; 26 k = 0; 27 for (i=0; i<mbs; i++) { 28 nz = ai[i+1] - ai[i]; 29 if (nz) { 30 rowlengths[k] += nz; /* no. of upper triangular blocks */ 31 if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */ 32 for (j=0; j<nz; j++) { /* no. of lower triangular blocks */ 33 rowlengths[(*aj)*bs]++; aj++; 34 } 35 } 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(((PetscObject)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_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 49 B->rmap->bs = A->rmap->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 - diagcnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz-diagcnt: %D\n",bi[mbs],2*a->nz - diagcnt); 66 67 /* set b->j and b->a */ 68 aj = a->j; av = a->a; 69 for (i=0; i<mbs; i++) { 70 nz = ai[i+1] - ai[i]; 71 /* diagonal block */ 72 if (nz && *aj == i) { 73 nz--; 74 for (j=0; j<bs; j++){ /* row i*bs+j */ 75 itmp = i*bs+j; 76 for (k=0; k<bs; k++){ /* col i*bs+k */ 77 *(bj + rowstart[itmp]) = (*aj)*bs+k; 78 *(bv + rowstart[itmp]) = *(av+k*bs+j); 79 rowstart[itmp]++; 80 } 81 } 82 aj++; av += bs2; 83 } 84 85 while (nz--){ 86 /* lower triangular blocks */ 87 for (j=0; j<bs; j++){ /* row (*aj)*bs+j */ 88 itmp = (*aj)*bs+j; 89 for (k=0; k<bs; k++){ /* col i*bs+k */ 90 *(bj + rowstart[itmp]) = i*bs+k; 91 *(bv + rowstart[itmp]) = *(av+j*bs+k); 92 rowstart[itmp]++; 93 } 94 } 95 /* upper triangular blocks */ 96 for (j=0; j<bs; j++){ /* row i*bs+j */ 97 itmp = i*bs+j; 98 for (k=0; k<bs; k++){ /* col (*aj)*bs+k */ 99 *(bj + rowstart[itmp]) = (*aj)*bs+k; 100 *(bv + rowstart[itmp]) = *(av+k*bs+j); 101 rowstart[itmp]++; 102 } 103 } 104 aj++; av += bs2; 105 } 106 } 107 ierr = PetscFree2(rowlengths,rowstart);CHKERRQ(ierr); 108 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 109 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 110 111 if (reuse == MAT_REUSE_MATRIX) { 112 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 113 } else { 114 *newmat = B; 115 } 116 PetscFunctionReturn(0); 117 } 118 EXTERN_C_END 119 120 EXTERN_C_BEGIN 121 #undef __FUNCT__ 122 #define __FUNCT__ "MatConvert_SeqAIJ_SeqSBAIJ" 123 PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,const MatType newtype,MatReuse reuse,Mat *newmat) { 124 Mat B; 125 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 126 Mat_SeqSBAIJ *b; 127 PetscErrorCode ierr; 128 PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths; 129 MatScalar *av,*bv; 130 131 PetscFunctionBegin; 132 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 133 134 ierr = PetscMalloc(m*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 135 for (i=0; i<m; i++) { 136 rowlengths[i] = ai[i+1] - a->diag[i]; 137 } 138 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 139 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 140 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 141 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,1,0,rowlengths);CHKERRQ(ierr); 142 143 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);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 if (A->hermitian){ 166 ierr = MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr); 167 } 168 169 if (reuse == MAT_REUSE_MATRIX) { 170 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 171 } else { 172 *newmat = B; 173 } 174 PetscFunctionReturn(0); 175 } 176 EXTERN_C_END 177 178 EXTERN_C_BEGIN 179 #undef __FUNCT__ 180 #define __FUNCT__ "MatConvert_SeqSBAIJ_SeqBAIJ" 181 PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 182 { 183 Mat B; 184 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 185 Mat_SeqBAIJ *b; 186 PetscErrorCode ierr; 187 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp; 188 PetscInt bs=A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row; 189 MatScalar *av,*bv; 190 191 PetscFunctionBegin; 192 /* compute browlengths of newmat */ 193 ierr = PetscMalloc2(mbs,PetscInt,&browlengths,mbs,PetscInt,&browstart);CHKERRQ(ierr); 194 for (i=0; i<mbs; i++) browlengths[i] = 0; 195 aj = a->j; 196 for (i=0; i<mbs; i++) { 197 nz = ai[i+1] - ai[i]; 198 aj++; /* skip diagonal */ 199 for (k=1; k<nz; k++) { /* no. of lower triangular blocks */ 200 browlengths[*aj]++; aj++; 201 } 202 browlengths[i] += nz; /* no. of upper triangular blocks */ 203 } 204 205 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 206 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 207 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 208 ierr = MatSeqBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 209 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);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_COMM_SELF,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 - transpose blocks of A */ 239 *(bj + browstart[*aj]) = i; /* block col index */ 240 itmp = bs2*browstart[*aj]; /* row index */ 241 for (col=0; col<bs; col++){ 242 k = col; 243 for (row=0; row<bs; row++){ 244 bv[itmp + col*bs+row] = av[k]; k+=bs; 245 } 246 } 247 browstart[*aj]++; 248 249 /* upper triangular blocks */ 250 *(bj + browstart[i]) = *aj; aj++; 251 itmp = bs2*browstart[i]; 252 for (k=0; k<bs2; k++){ 253 bv[itmp + k] = av[k]; 254 } 255 av += bs2; 256 browstart[i]++; 257 } 258 } 259 ierr = PetscFree2(browlengths,browstart);CHKERRQ(ierr); 260 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 261 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 262 263 if (reuse == MAT_REUSE_MATRIX) { 264 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 265 } else { 266 *newmat = B; 267 } 268 PetscFunctionReturn(0); 269 } 270 EXTERN_C_END 271 272 EXTERN_C_BEGIN 273 #undef __FUNCT__ 274 #define __FUNCT__ "MatConvert_SeqBAIJ_SeqSBAIJ" 275 PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 276 { 277 Mat B; 278 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 279 Mat_SeqSBAIJ *b; 280 PetscErrorCode ierr; 281 PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths; 282 PetscInt bs=A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd; 283 MatScalar *av,*bv; 284 PetscBool flg; 285 286 PetscFunctionBegin; 287 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 288 ierr = MatMissingDiagonal_SeqBAIJ(A,&flg,&dd);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */ 289 if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd); 290 291 ierr = PetscMalloc(mbs*sizeof(PetscInt),&browlengths);CHKERRQ(ierr); 292 for (i=0; i<mbs; i++) { 293 browlengths[i] = ai[i+1] - a->diag[i]; 294 } 295 296 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 297 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 298 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 299 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,browlengths);CHKERRQ(ierr); 300 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 301 302 b = (Mat_SeqSBAIJ*)(B->data); 303 bi = b->i; 304 bj = b->j; 305 bv = b->a; 306 307 bi[0] = 0; 308 for (i=0; i<mbs; i++) { 309 aj = a->j + a->diag[i]; 310 av = a->a + (a->diag[i])*bs2; 311 for (j=0; j<browlengths[i]; j++){ 312 *bj = *aj; bj++; aj++; 313 for (k=0; k<bs2; k++){ 314 *bv = *av; bv++; av++; 315 } 316 } 317 bi[i+1] = bi[i] + browlengths[i]; 318 b->ilen[i] = browlengths[i]; 319 } 320 ierr = PetscFree(browlengths);CHKERRQ(ierr); 321 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 322 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 323 324 if (reuse == MAT_REUSE_MATRIX) { 325 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 326 } else { 327 *newmat = B; 328 } 329 330 PetscFunctionReturn(0); 331 } 332 EXTERN_C_END 333