1 2 /* 3 Defines the basic matrix operations for the SBAIJ (compressed row) 4 matrix storage format. 5 */ 6 #include <../src/mat/impls/baij/seq/baij.h> /*I "petscmat.h" I*/ 7 #include <../src/mat/impls/sbaij/seq/sbaij.h> 8 #include <petscblaslapack.h> 9 10 #include <../src/mat/impls/sbaij/seq/relax.h> 11 #define USESHORT 12 #include <../src/mat/impls/sbaij/seq/relax.h> 13 14 extern PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat,PetscBool ); 15 16 /* 17 Checks for missing diagonals 18 */ 19 #undef __FUNCT__ 20 #define __FUNCT__ "MatMissingDiagonal_SeqSBAIJ" 21 PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A,PetscBool *missing,PetscInt *dd) 22 { 23 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 24 PetscErrorCode ierr; 25 PetscInt *diag,*jj = a->j,i; 26 27 PetscFunctionBegin; 28 ierr = MatMarkDiagonal_SeqSBAIJ(A);CHKERRQ(ierr); 29 diag = a->diag; 30 *missing = PETSC_FALSE; 31 for (i=0; i<a->mbs; i++) { 32 if (jj[diag[i]] != i) { 33 *missing = PETSC_TRUE; 34 if (dd) *dd = i; 35 break; 36 } 37 } 38 PetscFunctionReturn(0); 39 } 40 41 #undef __FUNCT__ 42 #define __FUNCT__ "MatMarkDiagonal_SeqSBAIJ" 43 PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A) 44 { 45 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 46 PetscErrorCode ierr; 47 PetscInt i; 48 49 PetscFunctionBegin; 50 if (!a->diag) { 51 ierr = PetscMalloc(a->mbs*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 52 ierr = PetscLogObjectMemory(A,a->mbs*sizeof(PetscInt));CHKERRQ(ierr); 53 a->free_diag = PETSC_TRUE; 54 } 55 for (i=0; i<a->mbs; i++) a->diag[i] = a->i[i]; 56 PetscFunctionReturn(0); 57 } 58 59 #undef __FUNCT__ 60 #define __FUNCT__ "MatGetRowIJ_SeqSBAIJ" 61 static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscBool *done) 62 { 63 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 64 PetscInt i,j,n = a->mbs,nz = a->i[n],bs = A->rmap->bs; 65 PetscErrorCode ierr; 66 67 PetscFunctionBegin; 68 *nn = n; 69 if (!ia) PetscFunctionReturn(0); 70 if (!blockcompressed) { 71 /* malloc & create the natural set of indices */ 72 ierr = PetscMalloc2((n+1)*bs,PetscInt,ia,nz*bs,PetscInt,ja);CHKERRQ(ierr); 73 for (i=0; i<n+1; i++) { 74 for (j=0; j<bs; j++) { 75 *ia[i*bs+j] = a->i[i]*bs+j+oshift; 76 } 77 } 78 for (i=0; i<nz; i++) { 79 for (j=0; j<bs; j++) { 80 *ja[i*bs+j] = a->j[i]*bs+j+oshift; 81 } 82 } 83 } else { /* blockcompressed */ 84 if (oshift == 1) { 85 /* temporarily add 1 to i and j indices */ 86 for (i=0; i<nz; i++) a->j[i]++; 87 for (i=0; i<n+1; i++) a->i[i]++; 88 } 89 *ia = a->i; *ja = a->j; 90 } 91 92 PetscFunctionReturn(0); 93 } 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "MatRestoreRowIJ_SeqSBAIJ" 97 static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscBool *done) 98 { 99 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 100 PetscInt i,n = a->mbs,nz = a->i[n]; 101 PetscErrorCode ierr; 102 103 PetscFunctionBegin; 104 if (!ia) PetscFunctionReturn(0); 105 106 if (!blockcompressed) { 107 ierr = PetscFree2(*ia,*ja);CHKERRQ(ierr); 108 } else if (oshift == 1) { /* blockcompressed */ 109 for (i=0; i<nz; i++) a->j[i]--; 110 for (i=0; i<n+1; i++) a->i[i]--; 111 } 112 113 PetscFunctionReturn(0); 114 } 115 116 #undef __FUNCT__ 117 #define __FUNCT__ "MatDestroy_SeqSBAIJ" 118 PetscErrorCode MatDestroy_SeqSBAIJ(Mat A) 119 { 120 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 121 PetscErrorCode ierr; 122 123 PetscFunctionBegin; 124 #if defined(PETSC_USE_LOG) 125 PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap->N,a->nz); 126 #endif 127 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 128 if (a->free_diag){ierr = PetscFree(a->diag);CHKERRQ(ierr);} 129 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 130 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 131 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 132 ierr = PetscFree(a->idiag);CHKERRQ(ierr); 133 ierr = PetscFree(a->inode.size);CHKERRQ(ierr); 134 ierr = PetscFree(a->diag);CHKERRQ(ierr); 135 if (a->free_imax_ilen) {ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr);} 136 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 137 ierr = PetscFree(a->sor_work);CHKERRQ(ierr); 138 ierr = PetscFree(a->solves_work);CHKERRQ(ierr); 139 ierr = PetscFree(a->mult_work);CHKERRQ(ierr); 140 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 141 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 142 if (a->free_jshort) {ierr = PetscFree(a->jshort);CHKERRQ(ierr);} 143 ierr = PetscFree(a->inew);CHKERRQ(ierr); 144 ierr = MatDestroy(&a->parent);CHKERRQ(ierr); 145 ierr = PetscFree(A->data);CHKERRQ(ierr); 146 147 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 148 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 149 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 150 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 151 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C","",PETSC_NULL);CHKERRQ(ierr); 152 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 153 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 154 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqsbstrm_C","",PETSC_NULL);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "MatSetOption_SeqSBAIJ" 160 PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscBool flg) 161 { 162 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 163 PetscErrorCode ierr; 164 165 PetscFunctionBegin; 166 switch (op) { 167 case MAT_ROW_ORIENTED: 168 a->roworiented = flg; 169 break; 170 case MAT_KEEP_NONZERO_PATTERN: 171 a->keepnonzeropattern = flg; 172 break; 173 case MAT_NEW_NONZERO_LOCATIONS: 174 a->nonew = (flg ? 0 : 1); 175 break; 176 case MAT_NEW_NONZERO_LOCATION_ERR: 177 a->nonew = (flg ? -1 : 0); 178 break; 179 case MAT_NEW_NONZERO_ALLOCATION_ERR: 180 a->nonew = (flg ? -2 : 0); 181 break; 182 case MAT_UNUSED_NONZERO_LOCATION_ERR: 183 a->nounused = (flg ? -1 : 0); 184 break; 185 case MAT_NEW_DIAGONALS: 186 case MAT_IGNORE_OFF_PROC_ENTRIES: 187 case MAT_USE_HASH_TABLE: 188 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 189 break; 190 case MAT_HERMITIAN: 191 if (!A->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatAssemblyEnd() first"); 192 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 193 A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian_ushort; 194 } else if (A->cmap->bs == 1) { 195 A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian; 196 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for Hermitian with block size greater than 1"); 197 break; 198 case MAT_SPD: 199 A->spd_set = PETSC_TRUE; 200 A->spd = flg; 201 if (flg) { 202 A->symmetric = PETSC_TRUE; 203 A->structurally_symmetric = PETSC_TRUE; 204 A->symmetric_set = PETSC_TRUE; 205 A->structurally_symmetric_set = PETSC_TRUE; 206 } 207 break; 208 case MAT_SYMMETRIC: 209 case MAT_STRUCTURALLY_SYMMETRIC: 210 case MAT_SYMMETRY_ETERNAL: 211 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix must be symmetric"); 212 ierr = PetscInfo1(A,"Option %s not relevent\n",MatOptions[op]);CHKERRQ(ierr); 213 break; 214 case MAT_IGNORE_LOWER_TRIANGULAR: 215 a->ignore_ltriangular = flg; 216 break; 217 case MAT_ERROR_LOWER_TRIANGULAR: 218 a->ignore_ltriangular = flg; 219 break; 220 case MAT_GETROW_UPPERTRIANGULAR: 221 a->getrow_utriangular = flg; 222 break; 223 default: 224 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 225 } 226 PetscFunctionReturn(0); 227 } 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "MatGetRow_SeqSBAIJ" 231 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **v) 232 { 233 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 234 PetscErrorCode ierr; 235 PetscInt itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2; 236 MatScalar *aa,*aa_i; 237 PetscScalar *v_i; 238 239 PetscFunctionBegin; 240 if (A && !a->getrow_utriangular) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatGetRow is not supported for SBAIJ matrix format. Getting the upper triangular part of row, run with -mat_getrow_uppertriangular, call MatSetOption(mat,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE) or MatGetRowUpperTriangular()"); 241 /* Get the upper triangular part of the row */ 242 bs = A->rmap->bs; 243 ai = a->i; 244 aj = a->j; 245 aa = a->a; 246 bs2 = a->bs2; 247 248 if (row < 0 || row >= A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Row %D out of range", row); 249 250 bn = row/bs; /* Block number */ 251 bp = row % bs; /* Block position */ 252 M = ai[bn+1] - ai[bn]; 253 *ncols = bs*M; 254 255 if (v) { 256 *v = 0; 257 if (*ncols) { 258 ierr = PetscMalloc((*ncols+row)*sizeof(PetscScalar),v);CHKERRQ(ierr); 259 for (i=0; i<M; i++) { /* for each block in the block row */ 260 v_i = *v + i*bs; 261 aa_i = aa + bs2*(ai[bn] + i); 262 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 263 } 264 } 265 } 266 267 if (cols) { 268 *cols = 0; 269 if (*ncols) { 270 ierr = PetscMalloc((*ncols+row)*sizeof(PetscInt),cols);CHKERRQ(ierr); 271 for (i=0; i<M; i++) { /* for each block in the block row */ 272 cols_i = *cols + i*bs; 273 itmp = bs*aj[ai[bn] + i]; 274 for (j=0; j<bs; j++) {cols_i[j] = itmp++;} 275 } 276 } 277 } 278 279 /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */ 280 /* this segment is currently removed, so only entries in the upper triangle are obtained */ 281 #ifdef column_search 282 v_i = *v + M*bs; 283 cols_i = *cols + M*bs; 284 for (i=0; i<bn; i++){ /* for each block row */ 285 M = ai[i+1] - ai[i]; 286 for (j=0; j<M; j++){ 287 itmp = aj[ai[i] + j]; /* block column value */ 288 if (itmp == bn){ 289 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 290 for (k=0; k<bs; k++) { 291 *cols_i++ = i*bs+k; 292 *v_i++ = aa_i[k]; 293 } 294 *ncols += bs; 295 break; 296 } 297 } 298 } 299 #endif 300 PetscFunctionReturn(0); 301 } 302 303 #undef __FUNCT__ 304 #define __FUNCT__ "MatRestoreRow_SeqSBAIJ" 305 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 306 { 307 PetscErrorCode ierr; 308 309 PetscFunctionBegin; 310 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);} 311 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);} 312 PetscFunctionReturn(0); 313 } 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "MatGetRowUpperTriangular_SeqSBAIJ" 317 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 318 { 319 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 320 321 PetscFunctionBegin; 322 a->getrow_utriangular = PETSC_TRUE; 323 PetscFunctionReturn(0); 324 } 325 #undef __FUNCT__ 326 #define __FUNCT__ "MatRestoreRowUpperTriangular_SeqSBAIJ" 327 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) 328 { 329 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 330 331 PetscFunctionBegin; 332 a->getrow_utriangular = PETSC_FALSE; 333 PetscFunctionReturn(0); 334 } 335 336 #undef __FUNCT__ 337 #define __FUNCT__ "MatTranspose_SeqSBAIJ" 338 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B) 339 { 340 PetscErrorCode ierr; 341 PetscFunctionBegin; 342 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 343 ierr = MatDuplicate(A,MAT_COPY_VALUES,B);CHKERRQ(ierr); 344 } 345 PetscFunctionReturn(0); 346 } 347 348 #undef __FUNCT__ 349 #define __FUNCT__ "MatView_SeqSBAIJ_ASCII" 350 static PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 351 { 352 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 353 PetscErrorCode ierr; 354 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2; 355 PetscViewerFormat format; 356 PetscInt *diag; 357 358 PetscFunctionBegin; 359 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 360 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 361 ierr = PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);CHKERRQ(ierr); 362 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 363 Mat aij; 364 if (A->factortype && bs>1){ 365 ierr = PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");CHKERRQ(ierr); 366 PetscFunctionReturn(0); 367 } 368 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);CHKERRQ(ierr); 369 ierr = MatView(aij,viewer);CHKERRQ(ierr); 370 ierr = MatDestroy(&aij);CHKERRQ(ierr); 371 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 372 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 373 for (i=0; i<a->mbs; i++) { 374 for (j=0; j<bs; j++) { 375 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 376 for (k=a->i[i]; k<a->i[i+1]; k++) { 377 for (l=0; l<bs; l++) { 378 #if defined(PETSC_USE_COMPLEX) 379 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 380 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 381 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 382 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 383 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 384 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 385 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 386 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 387 } 388 #else 389 if (a->a[bs2*k + l*bs + j] != 0.0) { 390 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 391 } 392 #endif 393 } 394 } 395 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 396 } 397 } 398 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 399 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 400 PetscFunctionReturn(0); 401 } else { 402 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 403 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 404 if (A->factortype){ /* for factored matrix */ 405 if (bs>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet"); 406 407 diag=a->diag; 408 for (i=0; i<a->mbs; i++) { /* for row block i */ 409 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 410 /* diagonal entry */ 411 #if defined(PETSC_USE_COMPLEX) 412 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) { 413 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",a->j[diag[i]], 414 PetscRealPart(1.0/a->a[diag[i]]),PetscImaginaryPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 415 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) { 416 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",a->j[diag[i]], 417 PetscRealPart(1.0/a->a[diag[i]]),-PetscImaginaryPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 418 } else { 419 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[diag[i]],PetscRealPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 420 } 421 #else 422 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[diag[i]],1.0/a->a[diag[i]]);CHKERRQ(ierr); 423 #endif 424 /* off-diagonal entries */ 425 for (k=a->i[i]; k<a->i[i+1]-1; k++) { 426 #if defined(PETSC_USE_COMPLEX) 427 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 428 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 429 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 430 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 431 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 432 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 433 } else { 434 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 435 } 436 #else 437 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[k],a->a[k]);CHKERRQ(ierr); 438 #endif 439 } 440 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 441 } 442 443 } else { /* for non-factored matrix */ 444 for (i=0; i<a->mbs; i++) { // for row block i 445 for (j=0; j<bs; j++) { // for row bs*i + j 446 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 447 for (k=a->i[i]; k<a->i[i+1]; k++) { // for column block 448 for (l=0; l<bs; l++) { // for column 449 #if defined(PETSC_USE_COMPLEX) 450 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 451 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 452 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 453 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 454 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 455 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 456 } else { 457 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 458 } 459 #else 460 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 461 #endif 462 } 463 } 464 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 465 } 466 } 467 } 468 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 469 } 470 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 471 PetscFunctionReturn(0); 472 } 473 474 #undef __FUNCT__ 475 #define __FUNCT__ "MatView_SeqSBAIJ_Draw_Zoom" 476 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 477 { 478 Mat A = (Mat) Aa; 479 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 480 PetscErrorCode ierr; 481 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 482 PetscMPIInt rank; 483 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 484 MatScalar *aa; 485 MPI_Comm comm; 486 PetscViewer viewer; 487 488 PetscFunctionBegin; 489 /* 490 This is nasty. If this is called from an originally parallel matrix 491 then all processes call this,but only the first has the matrix so the 492 rest should return immediately. 493 */ 494 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 495 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 496 if (rank) PetscFunctionReturn(0); 497 498 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 499 500 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 501 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 502 503 /* loop over matrix elements drawing boxes */ 504 color = PETSC_DRAW_BLUE; 505 for (i=0,row=0; i<mbs; i++,row+=bs) { 506 for (j=a->i[i]; j<a->i[i+1]; j++) { 507 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 508 x_l = a->j[j]*bs; x_r = x_l + 1.0; 509 aa = a->a + j*bs2; 510 for (k=0; k<bs; k++) { 511 for (l=0; l<bs; l++) { 512 if (PetscRealPart(*aa++) >= 0.) continue; 513 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 514 } 515 } 516 } 517 } 518 color = PETSC_DRAW_CYAN; 519 for (i=0,row=0; i<mbs; i++,row+=bs) { 520 for (j=a->i[i]; j<a->i[i+1]; j++) { 521 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 522 x_l = a->j[j]*bs; x_r = x_l + 1.0; 523 aa = a->a + j*bs2; 524 for (k=0; k<bs; k++) { 525 for (l=0; l<bs; l++) { 526 if (PetscRealPart(*aa++) != 0.) continue; 527 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 528 } 529 } 530 } 531 } 532 533 color = PETSC_DRAW_RED; 534 for (i=0,row=0; i<mbs; i++,row+=bs) { 535 for (j=a->i[i]; j<a->i[i+1]; j++) { 536 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 537 x_l = a->j[j]*bs; x_r = x_l + 1.0; 538 aa = a->a + j*bs2; 539 for (k=0; k<bs; k++) { 540 for (l=0; l<bs; l++) { 541 if (PetscRealPart(*aa++) <= 0.) continue; 542 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 543 } 544 } 545 } 546 } 547 PetscFunctionReturn(0); 548 } 549 550 #undef __FUNCT__ 551 #define __FUNCT__ "MatView_SeqSBAIJ_Draw" 552 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 553 { 554 PetscErrorCode ierr; 555 PetscReal xl,yl,xr,yr,w,h; 556 PetscDraw draw; 557 PetscBool isnull; 558 559 PetscFunctionBegin; 560 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 561 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 562 563 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 564 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 565 xr += w; yr += h; xl = -w; yl = -h; 566 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 567 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 568 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 569 PetscFunctionReturn(0); 570 } 571 572 #undef __FUNCT__ 573 #define __FUNCT__ "MatView_SeqSBAIJ" 574 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 575 { 576 PetscErrorCode ierr; 577 PetscBool iascii,isdraw; 578 FILE *file = 0; 579 580 PetscFunctionBegin; 581 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 582 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 583 if (iascii){ 584 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 585 } else if (isdraw) { 586 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 587 } else { 588 Mat B; 589 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 590 ierr = MatView(B,viewer);CHKERRQ(ierr); 591 ierr = MatDestroy(&B);CHKERRQ(ierr); 592 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 593 if (file) { 594 fprintf(file,"-matload_block_size %d\n",(int)A->rmap->bs); 595 } 596 } 597 PetscFunctionReturn(0); 598 } 599 600 601 #undef __FUNCT__ 602 #define __FUNCT__ "MatGetValues_SeqSBAIJ" 603 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 604 { 605 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 606 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 607 PetscInt *ai = a->i,*ailen = a->ilen; 608 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 609 MatScalar *ap,*aa = a->a; 610 611 PetscFunctionBegin; 612 for (k=0; k<m; k++) { /* loop over rows */ 613 row = im[k]; brow = row/bs; 614 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 615 if (row >= A->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 616 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 617 nrow = ailen[brow]; 618 for (l=0; l<n; l++) { /* loop over columns */ 619 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 620 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 621 col = in[l] ; 622 bcol = col/bs; 623 cidx = col%bs; 624 ridx = row%bs; 625 high = nrow; 626 low = 0; /* assume unsorted */ 627 while (high-low > 5) { 628 t = (low+high)/2; 629 if (rp[t] > bcol) high = t; 630 else low = t; 631 } 632 for (i=low; i<high; i++) { 633 if (rp[i] > bcol) break; 634 if (rp[i] == bcol) { 635 *v++ = ap[bs2*i+bs*cidx+ridx]; 636 goto finished; 637 } 638 } 639 *v++ = 0.0; 640 finished:; 641 } 642 } 643 PetscFunctionReturn(0); 644 } 645 646 647 #undef __FUNCT__ 648 #define __FUNCT__ "MatSetValuesBlocked_SeqSBAIJ" 649 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 650 { 651 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 652 PetscErrorCode ierr; 653 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 654 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen; 655 PetscInt *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 656 PetscBool roworiented=a->roworiented; 657 const PetscScalar *value = v; 658 MatScalar *ap,*aa = a->a,*bap; 659 660 PetscFunctionBegin; 661 if (roworiented) { 662 stepval = (n-1)*bs; 663 } else { 664 stepval = (m-1)*bs; 665 } 666 for (k=0; k<m; k++) { /* loop over added rows */ 667 row = im[k]; 668 if (row < 0) continue; 669 #if defined(PETSC_USE_DEBUG) 670 if (row >= a->mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1); 671 #endif 672 rp = aj + ai[row]; 673 ap = aa + bs2*ai[row]; 674 rmax = imax[row]; 675 nrow = ailen[row]; 676 low = 0; 677 high = nrow; 678 for (l=0; l<n; l++) { /* loop over added columns */ 679 if (in[l] < 0) continue; 680 col = in[l]; 681 #if defined(PETSC_USE_DEBUG) 682 if (col >= a->nbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1); 683 #endif 684 if (col < row) { 685 if (a->ignore_ltriangular) { 686 continue; /* ignore lower triangular block */ 687 } else { 688 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)"); 689 } 690 } 691 if (roworiented) { 692 value = v + k*(stepval+bs)*bs + l*bs; 693 } else { 694 value = v + l*(stepval+bs)*bs + k*bs; 695 } 696 if (col <= lastcol) low = 0; else high = nrow; 697 lastcol = col; 698 while (high-low > 7) { 699 t = (low+high)/2; 700 if (rp[t] > col) high = t; 701 else low = t; 702 } 703 for (i=low; i<high; i++) { 704 if (rp[i] > col) break; 705 if (rp[i] == col) { 706 bap = ap + bs2*i; 707 if (roworiented) { 708 if (is == ADD_VALUES) { 709 for (ii=0; ii<bs; ii++,value+=stepval) { 710 for (jj=ii; jj<bs2; jj+=bs) { 711 bap[jj] += *value++; 712 } 713 } 714 } else { 715 for (ii=0; ii<bs; ii++,value+=stepval) { 716 for (jj=ii; jj<bs2; jj+=bs) { 717 bap[jj] = *value++; 718 } 719 } 720 } 721 } else { 722 if (is == ADD_VALUES) { 723 for (ii=0; ii<bs; ii++,value+=stepval) { 724 for (jj=0; jj<bs; jj++) { 725 *bap++ += *value++; 726 } 727 } 728 } else { 729 for (ii=0; ii<bs; ii++,value+=stepval) { 730 for (jj=0; jj<bs; jj++) { 731 *bap++ = *value++; 732 } 733 } 734 } 735 } 736 goto noinsert2; 737 } 738 } 739 if (nonew == 1) goto noinsert2; 740 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 741 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 742 N = nrow++ - 1; high++; 743 /* shift up all the later entries in this row */ 744 for (ii=N; ii>=i; ii--) { 745 rp[ii+1] = rp[ii]; 746 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 747 } 748 if (N >= i) { 749 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 750 } 751 rp[i] = col; 752 bap = ap + bs2*i; 753 if (roworiented) { 754 for (ii=0; ii<bs; ii++,value+=stepval) { 755 for (jj=ii; jj<bs2; jj+=bs) { 756 bap[jj] = *value++; 757 } 758 } 759 } else { 760 for (ii=0; ii<bs; ii++,value+=stepval) { 761 for (jj=0; jj<bs; jj++) { 762 *bap++ = *value++; 763 } 764 } 765 } 766 noinsert2:; 767 low = i; 768 } 769 ailen[row] = nrow; 770 } 771 PetscFunctionReturn(0); 772 } 773 774 /* 775 This is not yet used 776 */ 777 #undef __FUNCT__ 778 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode" 779 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 780 { 781 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 782 PetscErrorCode ierr; 783 const PetscInt *ai = a->i, *aj = a->j,*cols; 784 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 785 PetscBool flag; 786 787 PetscFunctionBegin; 788 ierr = PetscMalloc(m*sizeof(PetscInt),&ns);CHKERRQ(ierr); 789 while (i < m){ 790 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 791 /* Limits the number of elements in a node to 'a->inode.limit' */ 792 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 793 nzy = ai[j+1] - ai[j]; 794 if (nzy != (nzx - j + i)) break; 795 ierr = PetscMemcmp(aj + ai[i] + j - i,aj + ai[j],nzy*sizeof(PetscInt),&flag);CHKERRQ(ierr); 796 if (!flag) break; 797 } 798 ns[node_count++] = blk_size; 799 i = j; 800 } 801 if (!a->inode.size && m && node_count > .9*m) { 802 ierr = PetscFree(ns);CHKERRQ(ierr); 803 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 804 } else { 805 a->inode.node_count = node_count; 806 ierr = PetscMalloc(node_count*sizeof(PetscInt),&a->inode.size);CHKERRQ(ierr); 807 ierr = PetscLogObjectMemory(A,node_count*sizeof(PetscInt));CHKERRQ(ierr); 808 ierr = PetscMemcpy(a->inode.size,ns,node_count*sizeof(PetscInt)); 809 ierr = PetscFree(ns);CHKERRQ(ierr); 810 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 811 812 /* count collections of adjacent columns in each inode */ 813 row = 0; 814 cnt = 0; 815 for (i=0; i<node_count; i++) { 816 cols = aj + ai[row] + a->inode.size[i]; 817 nz = ai[row+1] - ai[row] - a->inode.size[i]; 818 for (j=1; j<nz; j++) { 819 if (cols[j] != cols[j-1]+1) { 820 cnt++; 821 } 822 } 823 cnt++; 824 row += a->inode.size[i]; 825 } 826 ierr = PetscMalloc(2*cnt*sizeof(PetscInt),&counts);CHKERRQ(ierr); 827 cnt = 0; 828 row = 0; 829 for (i=0; i<node_count; i++) { 830 cols = aj + ai[row] + a->inode.size[i]; 831 CHKMEMQ; 832 counts[2*cnt] = cols[0]; 833 CHKMEMQ; 834 nz = ai[row+1] - ai[row] - a->inode.size[i]; 835 cnt2 = 1; 836 for (j=1; j<nz; j++) { 837 if (cols[j] != cols[j-1]+1) { 838 CHKMEMQ; 839 counts[2*(cnt++)+1] = cnt2; 840 counts[2*cnt] = cols[j]; 841 CHKMEMQ; 842 cnt2 = 1; 843 } else cnt2++; 844 } 845 CHKMEMQ; 846 counts[2*(cnt++)+1] = cnt2; 847 CHKMEMQ; 848 row += a->inode.size[i]; 849 } 850 ierr = PetscIntView(2*cnt,counts,0); 851 } 852 PetscFunctionReturn(0); 853 } 854 855 #undef __FUNCT__ 856 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 857 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 858 { 859 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 860 PetscErrorCode ierr; 861 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 862 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 863 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 864 MatScalar *aa = a->a,*ap; 865 866 PetscFunctionBegin; 867 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 868 869 if (m) rmax = ailen[0]; 870 for (i=1; i<mbs; i++) { 871 /* move each row back by the amount of empty slots (fshift) before it*/ 872 fshift += imax[i-1] - ailen[i-1]; 873 rmax = PetscMax(rmax,ailen[i]); 874 if (fshift) { 875 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 876 N = ailen[i]; 877 for (j=0; j<N; j++) { 878 ip[j-fshift] = ip[j]; 879 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 880 } 881 } 882 ai[i] = ai[i-1] + ailen[i-1]; 883 } 884 if (mbs) { 885 fshift += imax[mbs-1] - ailen[mbs-1]; 886 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 887 } 888 /* reset ilen and imax for each row */ 889 for (i=0; i<mbs; i++) { 890 ailen[i] = imax[i] = ai[i+1] - ai[i]; 891 } 892 a->nz = ai[mbs]; 893 894 /* diagonals may have moved, reset it */ 895 if (a->diag) { 896 ierr = PetscMemcpy(a->diag,ai,mbs*sizeof(PetscInt));CHKERRQ(ierr); 897 } 898 if (fshift && a->nounused == -1) { 899 SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D block size %D, %D unneeded", m, A->cmap->n, A->rmap->bs, fshift*bs2); 900 } 901 ierr = PetscInfo5(A,"Matrix size: %D X %D, block size %D; storage space: %D unneeded, %D used\n",m,A->rmap->N,A->rmap->bs,fshift*bs2,a->nz*bs2);CHKERRQ(ierr); 902 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 903 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 904 A->info.mallocs += a->reallocs; 905 a->reallocs = 0; 906 A->info.nz_unneeded = (PetscReal)fshift*bs2; 907 a->idiagvalid = PETSC_FALSE; 908 909 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 910 if (!a->jshort) { 911 ierr = PetscMalloc(a->i[A->rmap->n]*sizeof(unsigned short),&a->jshort);CHKERRQ(ierr); 912 ierr = PetscLogObjectMemory(A,a->i[A->rmap->n]*sizeof(unsigned short));CHKERRQ(ierr); 913 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 914 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 915 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 916 a->free_jshort = PETSC_TRUE; 917 } 918 } 919 PetscFunctionReturn(0); 920 } 921 922 /* 923 This function returns an array of flags which indicate the locations of contiguous 924 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 925 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 926 Assume: sizes should be long enough to hold all the values. 927 */ 928 #undef __FUNCT__ 929 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 930 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 931 { 932 PetscInt i,j,k,row; 933 PetscBool flg; 934 935 PetscFunctionBegin; 936 for (i=0,j=0; i<n; j++) { 937 row = idx[i]; 938 if (row%bs!=0) { /* Not the begining of a block */ 939 sizes[j] = 1; 940 i++; 941 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 942 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 943 i++; 944 } else { /* Begining of the block, so check if the complete block exists */ 945 flg = PETSC_TRUE; 946 for (k=1; k<bs; k++) { 947 if (row+k != idx[i+k]) { /* break in the block */ 948 flg = PETSC_FALSE; 949 break; 950 } 951 } 952 if (flg) { /* No break in the bs */ 953 sizes[j] = bs; 954 i+= bs; 955 } else { 956 sizes[j] = 1; 957 i++; 958 } 959 } 960 } 961 *bs_max = j; 962 PetscFunctionReturn(0); 963 } 964 965 966 /* Only add/insert a(i,j) with i<=j (blocks). 967 Any a(i,j) with i>j input by user is ingored. 968 */ 969 970 #undef __FUNCT__ 971 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 972 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 973 { 974 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 975 PetscErrorCode ierr; 976 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 977 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 978 PetscInt *aj=a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 979 PetscInt ridx,cidx,bs2=a->bs2; 980 MatScalar *ap,value,*aa=a->a,*bap; 981 982 PetscFunctionBegin; 983 if (v) PetscValidScalarPointer(v,6); 984 for (k=0; k<m; k++) { /* loop over added rows */ 985 row = im[k]; /* row number */ 986 brow = row/bs; /* block row number */ 987 if (row < 0) continue; 988 #if defined(PETSC_USE_DEBUG) 989 if (row >= A->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 990 #endif 991 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 992 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 993 rmax = imax[brow]; /* maximum space allocated for this row */ 994 nrow = ailen[brow]; /* actual length of this row */ 995 low = 0; 996 997 for (l=0; l<n; l++) { /* loop over added columns */ 998 if (in[l] < 0) continue; 999 #if defined(PETSC_USE_DEBUG) 1000 if (in[l] >= A->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->rmap->N-1); 1001 #endif 1002 col = in[l]; 1003 bcol = col/bs; /* block col number */ 1004 1005 if (brow > bcol) { 1006 if (a->ignore_ltriangular){ 1007 continue; /* ignore lower triangular values */ 1008 } else { 1009 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)"); 1010 } 1011 } 1012 1013 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 1014 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 1015 /* element value a(k,l) */ 1016 if (roworiented) { 1017 value = v[l + k*n]; 1018 } else { 1019 value = v[k + l*m]; 1020 } 1021 1022 /* move pointer bap to a(k,l) quickly and add/insert value */ 1023 if (col <= lastcol) low = 0; high = nrow; 1024 lastcol = col; 1025 while (high-low > 7) { 1026 t = (low+high)/2; 1027 if (rp[t] > bcol) high = t; 1028 else low = t; 1029 } 1030 for (i=low; i<high; i++) { 1031 if (rp[i] > bcol) break; 1032 if (rp[i] == bcol) { 1033 bap = ap + bs2*i + bs*cidx + ridx; 1034 if (is == ADD_VALUES) *bap += value; 1035 else *bap = value; 1036 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1037 if (brow == bcol && ridx < cidx){ 1038 bap = ap + bs2*i + bs*ridx + cidx; 1039 if (is == ADD_VALUES) *bap += value; 1040 else *bap = value; 1041 } 1042 goto noinsert1; 1043 } 1044 } 1045 1046 if (nonew == 1) goto noinsert1; 1047 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 1048 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 1049 1050 N = nrow++ - 1; high++; 1051 /* shift up all the later entries in this row */ 1052 for (ii=N; ii>=i; ii--) { 1053 rp[ii+1] = rp[ii]; 1054 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 1055 } 1056 if (N>=i) { 1057 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1058 } 1059 rp[i] = bcol; 1060 ap[bs2*i + bs*cidx + ridx] = value; 1061 noinsert1:; 1062 low = i; 1063 } 1064 } /* end of loop over added columns */ 1065 ailen[brow] = nrow; 1066 } /* end of loop over added rows */ 1067 PetscFunctionReturn(0); 1068 } 1069 1070 #undef __FUNCT__ 1071 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 1072 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1073 { 1074 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1075 Mat outA; 1076 PetscErrorCode ierr; 1077 PetscBool row_identity; 1078 1079 PetscFunctionBegin; 1080 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1081 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1082 if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1083 if (inA->rmap->bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix block size %D is not supported",inA->rmap->bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */ 1084 1085 outA = inA; 1086 inA->factortype = MAT_FACTOR_ICC; 1087 1088 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1089 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr); 1090 1091 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1092 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1093 a->row = row; 1094 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1095 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1096 a->col = row; 1097 1098 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1099 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1100 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1101 1102 if (!a->solve_work) { 1103 ierr = PetscMalloc((inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1104 ierr = PetscLogObjectMemory(inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1105 } 1106 1107 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1108 PetscFunctionReturn(0); 1109 } 1110 1111 EXTERN_C_BEGIN 1112 #undef __FUNCT__ 1113 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1114 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1115 { 1116 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1117 PetscInt i,nz,n; 1118 1119 PetscFunctionBegin; 1120 nz = baij->maxnz; 1121 n = mat->cmap->n; 1122 for (i=0; i<nz; i++) { 1123 baij->j[i] = indices[i]; 1124 } 1125 baij->nz = nz; 1126 for (i=0; i<n; i++) { 1127 baij->ilen[i] = baij->imax[i]; 1128 } 1129 PetscFunctionReturn(0); 1130 } 1131 EXTERN_C_END 1132 1133 #undef __FUNCT__ 1134 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1135 /*@ 1136 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1137 in the matrix. 1138 1139 Input Parameters: 1140 + mat - the SeqSBAIJ matrix 1141 - indices - the column indices 1142 1143 Level: advanced 1144 1145 Notes: 1146 This can be called if you have precomputed the nonzero structure of the 1147 matrix and want to provide it to the matrix object to improve the performance 1148 of the MatSetValues() operation. 1149 1150 You MUST have set the correct numbers of nonzeros per row in the call to 1151 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1152 1153 MUST be called before any calls to MatSetValues() 1154 1155 .seealso: MatCreateSeqSBAIJ 1156 @*/ 1157 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1158 { 1159 PetscErrorCode ierr; 1160 1161 PetscFunctionBegin; 1162 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1163 PetscValidPointer(indices,2); 1164 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 1165 PetscFunctionReturn(0); 1166 } 1167 1168 #undef __FUNCT__ 1169 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1170 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1171 { 1172 PetscErrorCode ierr; 1173 1174 PetscFunctionBegin; 1175 /* If the two matrices have the same copy implementation, use fast copy. */ 1176 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1177 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1178 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1179 1180 if (a->i[A->rmap->N] != b->i[B->rmap->N]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1181 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1182 } else { 1183 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1184 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1185 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1186 } 1187 PetscFunctionReturn(0); 1188 } 1189 1190 #undef __FUNCT__ 1191 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1192 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1193 { 1194 PetscErrorCode ierr; 1195 1196 PetscFunctionBegin; 1197 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,-PetscMax(A->rmap->bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1198 PetscFunctionReturn(0); 1199 } 1200 1201 #undef __FUNCT__ 1202 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1203 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1204 { 1205 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1206 PetscFunctionBegin; 1207 *array = a->a; 1208 PetscFunctionReturn(0); 1209 } 1210 1211 #undef __FUNCT__ 1212 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1213 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1214 { 1215 PetscFunctionBegin; 1216 PetscFunctionReturn(0); 1217 } 1218 1219 #undef __FUNCT__ 1220 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1221 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1222 { 1223 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1224 PetscErrorCode ierr; 1225 PetscInt i,bs=Y->rmap->bs,bs2=bs*bs,j; 1226 PetscBLASInt one = 1; 1227 1228 PetscFunctionBegin; 1229 if (str == SAME_NONZERO_PATTERN) { 1230 PetscScalar alpha = a; 1231 PetscBLASInt bnz = PetscBLASIntCast(x->nz*bs2); 1232 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1233 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1234 if (y->xtoy && y->XtoY != X) { 1235 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1236 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 1237 } 1238 if (!y->xtoy) { /* get xtoy */ 1239 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1240 y->XtoY = X; 1241 } 1242 for (i=0; i<x->nz; i++) { 1243 j = 0; 1244 while (j < bs2){ 1245 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1246 j++; 1247 } 1248 } 1249 ierr = PetscInfo3(Y,"ratio of nnz_s(X)/nnz_s(Y): %D/%D = %G\n",bs2*x->nz,bs2*y->nz,(PetscReal)(bs2*x->nz)/(bs2*y->nz));CHKERRQ(ierr); 1250 } else { 1251 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1252 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1253 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1254 } 1255 PetscFunctionReturn(0); 1256 } 1257 1258 #undef __FUNCT__ 1259 #define __FUNCT__ "MatSetBlockSize_SeqSBAIJ" 1260 PetscErrorCode MatSetBlockSize_SeqSBAIJ(Mat A,PetscInt bs) 1261 { 1262 PetscInt rbs,cbs; 1263 PetscErrorCode ierr; 1264 1265 PetscFunctionBegin; 1266 ierr = PetscLayoutGetBlockSize(A->rmap,&rbs);CHKERRQ(ierr); 1267 ierr = PetscLayoutGetBlockSize(A->cmap,&cbs);CHKERRQ(ierr); 1268 if (rbs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Attempt to set block size %d with SBAIJ %d",bs,rbs); 1269 if (cbs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Attempt to set block size %d with SBAIJ %d",bs,cbs); 1270 PetscFunctionReturn(0); 1271 } 1272 1273 #undef __FUNCT__ 1274 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1275 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1276 { 1277 PetscFunctionBegin; 1278 *flg = PETSC_TRUE; 1279 PetscFunctionReturn(0); 1280 } 1281 1282 #undef __FUNCT__ 1283 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1284 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1285 { 1286 PetscFunctionBegin; 1287 *flg = PETSC_TRUE; 1288 PetscFunctionReturn(0); 1289 } 1290 1291 #undef __FUNCT__ 1292 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1293 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1294 { 1295 PetscFunctionBegin; 1296 *flg = PETSC_FALSE; 1297 PetscFunctionReturn(0); 1298 } 1299 1300 #undef __FUNCT__ 1301 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1302 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1303 { 1304 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1305 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1306 MatScalar *aa = a->a; 1307 1308 PetscFunctionBegin; 1309 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1310 PetscFunctionReturn(0); 1311 } 1312 1313 #undef __FUNCT__ 1314 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1315 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1316 { 1317 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1318 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1319 MatScalar *aa = a->a; 1320 1321 PetscFunctionBegin; 1322 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1323 PetscFunctionReturn(0); 1324 } 1325 1326 #undef __FUNCT__ 1327 #define __FUNCT__ "MatZeroRowsColumns_SeqSBAIJ" 1328 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1329 { 1330 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1331 PetscErrorCode ierr; 1332 PetscInt i,j,k,count; 1333 PetscInt bs=A->rmap->bs,bs2=baij->bs2,row,col; 1334 PetscScalar zero = 0.0; 1335 MatScalar *aa; 1336 const PetscScalar *xx; 1337 PetscScalar *bb; 1338 PetscBool *zeroed,vecs = PETSC_FALSE; 1339 1340 PetscFunctionBegin; 1341 /* fix right hand side if needed */ 1342 if (x && b) { 1343 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1344 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1345 vecs = PETSC_TRUE; 1346 } 1347 A->same_nonzero = PETSC_TRUE; 1348 1349 /* zero the columns */ 1350 ierr = PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);CHKERRQ(ierr); 1351 ierr = PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));CHKERRQ(ierr); 1352 for (i=0; i<is_n; i++) { 1353 if (is_idx[i] < 0 || is_idx[i] >= A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range",is_idx[i]); 1354 zeroed[is_idx[i]] = PETSC_TRUE; 1355 } 1356 if (vecs) { 1357 for (i=0; i<A->rmap->N; i++) { 1358 row = i/bs; 1359 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1360 for (k=0; k<bs; k++) { 1361 col = bs*baij->j[j] + k; 1362 if (col <= i) continue; 1363 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1364 if (!zeroed[i] && zeroed[col]) { 1365 bb[i] -= aa[0]*xx[col]; 1366 } 1367 if (zeroed[i] && !zeroed[col]) { 1368 bb[col] -= aa[0]*xx[i]; 1369 } 1370 } 1371 } 1372 } 1373 for (i=0; i<is_n; i++) { 1374 bb[is_idx[i]] = diag*xx[is_idx[i]]; 1375 } 1376 } 1377 1378 for (i=0; i<A->rmap->N; i++) { 1379 if (!zeroed[i]) { 1380 row = i/bs; 1381 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1382 for (k=0; k<bs; k++) { 1383 col = bs*baij->j[j] + k; 1384 if (zeroed[col]) { 1385 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1386 aa[0] = 0.0; 1387 } 1388 } 1389 } 1390 } 1391 } 1392 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1393 if (vecs) { 1394 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1395 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1396 } 1397 1398 /* zero the rows */ 1399 for (i=0; i<is_n; i++) { 1400 row = is_idx[i]; 1401 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1402 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1403 for (k=0; k<count; k++) { 1404 aa[0] = zero; 1405 aa += bs; 1406 } 1407 if (diag != 0.0) { 1408 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1409 } 1410 } 1411 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1412 PetscFunctionReturn(0); 1413 } 1414 1415 /* -------------------------------------------------------------------*/ 1416 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1417 MatGetRow_SeqSBAIJ, 1418 MatRestoreRow_SeqSBAIJ, 1419 MatMult_SeqSBAIJ_N, 1420 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1421 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1422 MatMultAdd_SeqSBAIJ_N, 1423 0, 1424 0, 1425 0, 1426 /*10*/ 0, 1427 0, 1428 MatCholeskyFactor_SeqSBAIJ, 1429 MatSOR_SeqSBAIJ, 1430 MatTranspose_SeqSBAIJ, 1431 /*15*/ MatGetInfo_SeqSBAIJ, 1432 MatEqual_SeqSBAIJ, 1433 MatGetDiagonal_SeqSBAIJ, 1434 MatDiagonalScale_SeqSBAIJ, 1435 MatNorm_SeqSBAIJ, 1436 /*20*/ 0, 1437 MatAssemblyEnd_SeqSBAIJ, 1438 MatSetOption_SeqSBAIJ, 1439 MatZeroEntries_SeqSBAIJ, 1440 /*24*/ 0, 1441 0, 1442 0, 1443 0, 1444 0, 1445 /*29*/ MatSetUpPreallocation_SeqSBAIJ, 1446 0, 1447 0, 1448 MatGetArray_SeqSBAIJ, 1449 MatRestoreArray_SeqSBAIJ, 1450 /*34*/ MatDuplicate_SeqSBAIJ, 1451 0, 1452 0, 1453 0, 1454 MatICCFactor_SeqSBAIJ, 1455 /*39*/ MatAXPY_SeqSBAIJ, 1456 MatGetSubMatrices_SeqSBAIJ, 1457 MatIncreaseOverlap_SeqSBAIJ, 1458 MatGetValues_SeqSBAIJ, 1459 MatCopy_SeqSBAIJ, 1460 /*44*/ 0, 1461 MatScale_SeqSBAIJ, 1462 0, 1463 0, 1464 MatZeroRowsColumns_SeqSBAIJ, 1465 /*49*/ MatSetBlockSize_SeqSBAIJ, 1466 MatGetRowIJ_SeqSBAIJ, 1467 MatRestoreRowIJ_SeqSBAIJ, 1468 0, 1469 0, 1470 /*54*/ 0, 1471 0, 1472 0, 1473 0, 1474 MatSetValuesBlocked_SeqSBAIJ, 1475 /*59*/ MatGetSubMatrix_SeqSBAIJ, 1476 0, 1477 0, 1478 0, 1479 0, 1480 /*64*/ 0, 1481 0, 1482 0, 1483 0, 1484 0, 1485 /*69*/ MatGetRowMaxAbs_SeqSBAIJ, 1486 0, 1487 0, 1488 0, 1489 0, 1490 /*74*/ 0, 1491 0, 1492 0, 1493 0, 1494 0, 1495 /*79*/ 0, 1496 0, 1497 0, 1498 MatGetInertia_SeqSBAIJ, 1499 MatLoad_SeqSBAIJ, 1500 /*84*/ MatIsSymmetric_SeqSBAIJ, 1501 MatIsHermitian_SeqSBAIJ, 1502 MatIsStructurallySymmetric_SeqSBAIJ, 1503 0, 1504 0, 1505 /*89*/ 0, 1506 0, 1507 0, 1508 0, 1509 0, 1510 /*94*/ 0, 1511 0, 1512 0, 1513 0, 1514 0, 1515 /*99*/ 0, 1516 0, 1517 0, 1518 0, 1519 0, 1520 /*104*/0, 1521 MatRealPart_SeqSBAIJ, 1522 MatImaginaryPart_SeqSBAIJ, 1523 MatGetRowUpperTriangular_SeqSBAIJ, 1524 MatRestoreRowUpperTriangular_SeqSBAIJ, 1525 /*109*/0, 1526 0, 1527 0, 1528 0, 1529 MatMissingDiagonal_SeqSBAIJ, 1530 /*114*/0, 1531 0, 1532 0, 1533 0, 1534 0, 1535 /*119*/0, 1536 0, 1537 0, 1538 0 1539 }; 1540 1541 EXTERN_C_BEGIN 1542 #undef __FUNCT__ 1543 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1544 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1545 { 1546 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1547 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1548 PetscErrorCode ierr; 1549 1550 PetscFunctionBegin; 1551 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1552 1553 /* allocate space for values if not already there */ 1554 if (!aij->saved_values) { 1555 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1556 } 1557 1558 /* copy values over */ 1559 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1560 PetscFunctionReturn(0); 1561 } 1562 EXTERN_C_END 1563 1564 EXTERN_C_BEGIN 1565 #undef __FUNCT__ 1566 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1567 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1568 { 1569 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1570 PetscErrorCode ierr; 1571 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1572 1573 PetscFunctionBegin; 1574 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1575 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1576 1577 /* copy values over */ 1578 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1579 PetscFunctionReturn(0); 1580 } 1581 EXTERN_C_END 1582 1583 EXTERN_C_BEGIN 1584 #undef __FUNCT__ 1585 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1586 PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1587 { 1588 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1589 PetscErrorCode ierr; 1590 PetscInt i,mbs,bs2, newbs = PetscAbs(bs); 1591 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE; 1592 1593 PetscFunctionBegin; 1594 B->preallocated = PETSC_TRUE; 1595 if (bs < 0) { 1596 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for MPISBAIJ matrix","Mat");CHKERRQ(ierr); 1597 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSeqSBAIJSetPreallocation",newbs,&newbs,PETSC_NULL);CHKERRQ(ierr); 1598 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1599 bs = PetscAbs(bs); 1600 } 1601 if (nnz && newbs != bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot change blocksize from command line if setting nnz"); 1602 bs = newbs; 1603 1604 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1605 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1606 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1607 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1608 1609 mbs = B->rmap->N/bs; 1610 bs2 = bs*bs; 1611 1612 if (mbs*bs != B->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1613 1614 if (nz == MAT_SKIP_ALLOCATION) { 1615 skipallocation = PETSC_TRUE; 1616 nz = 0; 1617 } 1618 1619 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1620 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1621 if (nnz) { 1622 for (i=0; i<mbs; i++) { 1623 if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1624 if (nnz[i] > mbs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %D value %D rowlength %D",i,nnz[i],mbs); 1625 } 1626 } 1627 1628 B->ops->mult = MatMult_SeqSBAIJ_N; 1629 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1630 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1631 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1632 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1633 if (!flg) { 1634 switch (bs) { 1635 case 1: 1636 B->ops->mult = MatMult_SeqSBAIJ_1; 1637 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1638 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1639 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1640 break; 1641 case 2: 1642 B->ops->mult = MatMult_SeqSBAIJ_2; 1643 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1644 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1645 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1646 break; 1647 case 3: 1648 B->ops->mult = MatMult_SeqSBAIJ_3; 1649 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1650 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1651 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1652 break; 1653 case 4: 1654 B->ops->mult = MatMult_SeqSBAIJ_4; 1655 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1656 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1657 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1658 break; 1659 case 5: 1660 B->ops->mult = MatMult_SeqSBAIJ_5; 1661 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1662 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1663 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1664 break; 1665 case 6: 1666 B->ops->mult = MatMult_SeqSBAIJ_6; 1667 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1668 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1669 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1670 break; 1671 case 7: 1672 B->ops->mult = MatMult_SeqSBAIJ_7; 1673 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1674 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1675 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1676 break; 1677 } 1678 } 1679 1680 b->mbs = mbs; 1681 b->nbs = mbs; 1682 if (!skipallocation) { 1683 if (!b->imax) { 1684 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1685 b->free_imax_ilen = PETSC_TRUE; 1686 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1687 } 1688 if (!nnz) { 1689 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1690 else if (nz <= 0) nz = 1; 1691 for (i=0; i<mbs; i++) { 1692 b->imax[i] = nz; 1693 } 1694 nz = nz*mbs; /* total nz */ 1695 } else { 1696 nz = 0; 1697 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1698 } 1699 /* b->ilen will count nonzeros in each block row so far. */ 1700 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1701 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1702 1703 /* allocate the matrix space */ 1704 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1705 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1706 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1707 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1708 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1709 b->singlemalloc = PETSC_TRUE; 1710 1711 /* pointer to beginning of each row */ 1712 b->i[0] = 0; 1713 for (i=1; i<mbs+1; i++) { 1714 b->i[i] = b->i[i-1] + b->imax[i-1]; 1715 } 1716 b->free_a = PETSC_TRUE; 1717 b->free_ij = PETSC_TRUE; 1718 } else { 1719 b->free_a = PETSC_FALSE; 1720 b->free_ij = PETSC_FALSE; 1721 } 1722 1723 B->rmap->bs = bs; 1724 b->bs2 = bs2; 1725 b->nz = 0; 1726 b->maxnz = nz; 1727 1728 b->inew = 0; 1729 b->jnew = 0; 1730 b->anew = 0; 1731 b->a2anew = 0; 1732 b->permute = PETSC_FALSE; 1733 PetscFunctionReturn(0); 1734 } 1735 EXTERN_C_END 1736 1737 /* 1738 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1739 */ 1740 #undef __FUNCT__ 1741 #define __FUNCT__ "MatSeqSBAIJSetNumericFactorization_inplace" 1742 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1743 { 1744 PetscErrorCode ierr; 1745 PetscBool flg = PETSC_FALSE; 1746 PetscInt bs = B->rmap->bs; 1747 1748 PetscFunctionBegin; 1749 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1750 if (flg) bs = 8; 1751 1752 if (!natural) { 1753 switch (bs) { 1754 case 1: 1755 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1756 break; 1757 case 2: 1758 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1759 break; 1760 case 3: 1761 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1762 break; 1763 case 4: 1764 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1765 break; 1766 case 5: 1767 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1768 break; 1769 case 6: 1770 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1771 break; 1772 case 7: 1773 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1774 break; 1775 default: 1776 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1777 break; 1778 } 1779 } else { 1780 switch (bs) { 1781 case 1: 1782 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1783 break; 1784 case 2: 1785 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1786 break; 1787 case 3: 1788 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1789 break; 1790 case 4: 1791 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1792 break; 1793 case 5: 1794 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1795 break; 1796 case 6: 1797 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1798 break; 1799 case 7: 1800 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1801 break; 1802 default: 1803 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1804 break; 1805 } 1806 } 1807 PetscFunctionReturn(0); 1808 } 1809 1810 EXTERN_C_BEGIN 1811 extern PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1812 extern PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1813 EXTERN_C_END 1814 1815 1816 EXTERN_C_BEGIN 1817 #undef __FUNCT__ 1818 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1819 PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1820 { 1821 PetscInt n = A->rmap->n; 1822 PetscErrorCode ierr; 1823 1824 PetscFunctionBegin; 1825 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1826 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1827 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1828 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1829 ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1830 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1831 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1832 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1833 (*B)->factortype = ftype; 1834 PetscFunctionReturn(0); 1835 } 1836 EXTERN_C_END 1837 1838 EXTERN_C_BEGIN 1839 #undef __FUNCT__ 1840 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1841 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscBool *flg) 1842 { 1843 PetscFunctionBegin; 1844 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1845 *flg = PETSC_TRUE; 1846 } else { 1847 *flg = PETSC_FALSE; 1848 } 1849 PetscFunctionReturn(0); 1850 } 1851 EXTERN_C_END 1852 1853 EXTERN_C_BEGIN 1854 #if defined(PETSC_HAVE_MUMPS) 1855 extern PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 1856 #endif 1857 #if defined(PETSC_HAVE_SPOOLES) 1858 extern PetscErrorCode MatGetFactor_seqsbaij_spooles(Mat,MatFactorType,Mat*); 1859 #endif 1860 #if defined(PETSC_HAVE_PASTIX) 1861 extern PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1862 #endif 1863 #if defined(PETSC_HAVE_CHOLMOD) 1864 extern PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*); 1865 #endif 1866 extern PetscErrorCode MatGetFactor_seqsbaij_sbstrm(Mat,MatFactorType,Mat*); 1867 EXTERN_C_END 1868 1869 /*MC 1870 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1871 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1872 1873 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1874 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1875 1876 Options Database Keys: 1877 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1878 1879 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1880 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1881 the options database -mat_ignore_lower_triangular false it will generate an error if you try to set a value in the lower triangular portion. 1882 1883 1884 Level: beginner 1885 1886 .seealso: MatCreateSeqSBAIJ 1887 M*/ 1888 1889 EXTERN_C_BEGIN 1890 extern PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1891 EXTERN_C_END 1892 1893 1894 EXTERN_C_BEGIN 1895 #undef __FUNCT__ 1896 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1897 PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1898 { 1899 Mat_SeqSBAIJ *b; 1900 PetscErrorCode ierr; 1901 PetscMPIInt size; 1902 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1903 1904 PetscFunctionBegin; 1905 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1906 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1907 1908 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1909 B->data = (void*)b; 1910 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1911 B->ops->destroy = MatDestroy_SeqSBAIJ; 1912 B->ops->view = MatView_SeqSBAIJ; 1913 b->row = 0; 1914 b->icol = 0; 1915 b->reallocs = 0; 1916 b->saved_values = 0; 1917 b->inode.limit = 5; 1918 b->inode.max_limit = 5; 1919 1920 b->roworiented = PETSC_TRUE; 1921 b->nonew = 0; 1922 b->diag = 0; 1923 b->solve_work = 0; 1924 b->mult_work = 0; 1925 B->spptr = 0; 1926 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1927 b->keepnonzeropattern = PETSC_FALSE; 1928 b->xtoy = 0; 1929 b->XtoY = 0; 1930 1931 b->inew = 0; 1932 b->jnew = 0; 1933 b->anew = 0; 1934 b->a2anew = 0; 1935 b->permute = PETSC_FALSE; 1936 1937 b->ignore_ltriangular = PETSC_TRUE; 1938 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,PETSC_NULL);CHKERRQ(ierr); 1939 1940 b->getrow_utriangular = PETSC_FALSE; 1941 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,PETSC_NULL);CHKERRQ(ierr); 1942 1943 #if defined(PETSC_HAVE_PASTIX) 1944 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 1945 "MatGetFactor_seqsbaij_pastix", 1946 MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1947 #endif 1948 #if defined(PETSC_HAVE_SPOOLES) 1949 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 1950 "MatGetFactor_seqsbaij_spooles", 1951 MatGetFactor_seqsbaij_spooles);CHKERRQ(ierr); 1952 #endif 1953 #if defined(PETSC_HAVE_MUMPS) 1954 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 1955 "MatGetFactor_sbaij_mumps", 1956 MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 1957 #endif 1958 #if defined(PETSC_HAVE_CHOLMOD) 1959 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 1960 "MatGetFactor_seqsbaij_cholmod", 1961 MatGetFactor_seqsbaij_cholmod);CHKERRQ(ierr); 1962 #endif 1963 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 1964 "MatGetFactorAvailable_seqsbaij_petsc", 1965 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1966 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 1967 "MatGetFactor_seqsbaij_petsc", 1968 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1969 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_sbstrm_C", 1970 "MatGetFactor_seqsbaij_sbstrm", 1971 MatGetFactor_seqsbaij_sbstrm);CHKERRQ(ierr); 1972 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1973 "MatStoreValues_SeqSBAIJ", 1974 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1975 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1976 "MatRetrieveValues_SeqSBAIJ", 1977 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1978 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1979 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1980 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1981 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1982 "MatConvert_SeqSBAIJ_SeqAIJ", 1983 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1984 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1985 "MatConvert_SeqSBAIJ_SeqBAIJ", 1986 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1987 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1988 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1989 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1990 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqsbstrm_C", 1991 "MatConvert_SeqSBAIJ_SeqSBSTRM", 1992 MatConvert_SeqSBAIJ_SeqSBSTRM);CHKERRQ(ierr); 1993 1994 B->symmetric = PETSC_TRUE; 1995 B->structurally_symmetric = PETSC_TRUE; 1996 B->symmetric_set = PETSC_TRUE; 1997 B->structurally_symmetric_set = PETSC_TRUE; 1998 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1999 2000 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 2001 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",PETSC_NULL,no_unroll,&no_unroll,PETSC_NULL);CHKERRQ(ierr); 2002 if (no_unroll) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr);} 2003 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",PETSC_NULL,no_inode,&no_inode,PETSC_NULL);CHKERRQ(ierr); 2004 if (no_inode) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr);} 2005 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",PETSC_NULL,b->inode.limit,&b->inode.limit,PETSC_NULL);CHKERRQ(ierr); 2006 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2007 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 2008 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 2009 2010 PetscFunctionReturn(0); 2011 } 2012 EXTERN_C_END 2013 2014 #undef __FUNCT__ 2015 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 2016 /*@C 2017 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2018 compressed row) format. For good matrix assembly performance the 2019 user should preallocate the matrix storage by setting the parameter nz 2020 (or the array nnz). By setting these parameters accurately, performance 2021 during matrix assembly can be increased by more than a factor of 50. 2022 2023 Collective on Mat 2024 2025 Input Parameters: 2026 + A - the symmetric matrix 2027 . bs - size of block 2028 . nz - number of block nonzeros per block row (same for all rows) 2029 - nnz - array containing the number of block nonzeros in the upper triangular plus 2030 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 2031 2032 Options Database Keys: 2033 . -mat_no_unroll - uses code that does not unroll the loops in the 2034 block calculations (much slower) 2035 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2036 2037 Level: intermediate 2038 2039 Notes: 2040 Specify the preallocated storage with either nz or nnz (not both). 2041 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2042 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2043 2044 You can call MatGetInfo() to get information on how effective the preallocation was; 2045 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2046 You can also run with the option -info and look for messages with the string 2047 malloc in them to see if additional memory allocation was needed. 2048 2049 If the nnz parameter is given then the nz parameter is ignored 2050 2051 2052 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 2053 @*/ 2054 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2055 { 2056 PetscErrorCode ierr; 2057 2058 PetscFunctionBegin; 2059 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2060 PetscFunctionReturn(0); 2061 } 2062 2063 #undef __FUNCT__ 2064 #define __FUNCT__ "MatCreateSeqSBAIJ" 2065 /*@C 2066 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2067 compressed row) format. For good matrix assembly performance the 2068 user should preallocate the matrix storage by setting the parameter nz 2069 (or the array nnz). By setting these parameters accurately, performance 2070 during matrix assembly can be increased by more than a factor of 50. 2071 2072 Collective on MPI_Comm 2073 2074 Input Parameters: 2075 + comm - MPI communicator, set to PETSC_COMM_SELF 2076 . bs - size of block 2077 . m - number of rows, or number of columns 2078 . nz - number of block nonzeros per block row (same for all rows) 2079 - nnz - array containing the number of block nonzeros in the upper triangular plus 2080 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 2081 2082 Output Parameter: 2083 . A - the symmetric matrix 2084 2085 Options Database Keys: 2086 . -mat_no_unroll - uses code that does not unroll the loops in the 2087 block calculations (much slower) 2088 . -mat_block_size - size of the blocks to use 2089 2090 Level: intermediate 2091 2092 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2093 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2094 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2095 2096 Notes: 2097 The number of rows and columns must be divisible by blocksize. 2098 This matrix type does not support complex Hermitian operation. 2099 2100 Specify the preallocated storage with either nz or nnz (not both). 2101 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2102 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2103 2104 If the nnz parameter is given then the nz parameter is ignored 2105 2106 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 2107 @*/ 2108 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2109 { 2110 PetscErrorCode ierr; 2111 2112 PetscFunctionBegin; 2113 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2114 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2115 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2116 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2117 PetscFunctionReturn(0); 2118 } 2119 2120 #undef __FUNCT__ 2121 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2122 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2123 { 2124 Mat C; 2125 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2126 PetscErrorCode ierr; 2127 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2128 2129 PetscFunctionBegin; 2130 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2131 2132 *B = 0; 2133 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 2134 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2135 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2136 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2137 c = (Mat_SeqSBAIJ*)C->data; 2138 2139 C->preallocated = PETSC_TRUE; 2140 C->factortype = A->factortype; 2141 c->row = 0; 2142 c->icol = 0; 2143 c->saved_values = 0; 2144 c->keepnonzeropattern = a->keepnonzeropattern; 2145 C->assembled = PETSC_TRUE; 2146 2147 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2148 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2149 c->bs2 = a->bs2; 2150 c->mbs = a->mbs; 2151 c->nbs = a->nbs; 2152 2153 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2154 c->imax = a->imax; 2155 c->ilen = a->ilen; 2156 c->free_imax_ilen = PETSC_FALSE; 2157 } else { 2158 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 2159 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2160 for (i=0; i<mbs; i++) { 2161 c->imax[i] = a->imax[i]; 2162 c->ilen[i] = a->ilen[i]; 2163 } 2164 c->free_imax_ilen = PETSC_TRUE; 2165 } 2166 2167 /* allocate the matrix space */ 2168 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2169 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2170 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2171 c->singlemalloc = PETSC_FALSE; 2172 c->free_ij = PETSC_FALSE; 2173 c->parent = A; 2174 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2175 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2176 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2177 } else { 2178 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2179 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2180 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2181 c->singlemalloc = PETSC_TRUE; 2182 c->free_ij = PETSC_TRUE; 2183 } 2184 if (mbs > 0) { 2185 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2186 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2187 } 2188 if (cpvalues == MAT_COPY_VALUES) { 2189 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2190 } else { 2191 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2192 } 2193 if (a->jshort) { 2194 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2195 c->jshort = a->jshort; 2196 c->free_jshort = PETSC_FALSE; 2197 } else { 2198 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2199 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2200 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2201 c->free_jshort = PETSC_TRUE; 2202 } 2203 } 2204 } 2205 2206 c->roworiented = a->roworiented; 2207 c->nonew = a->nonew; 2208 2209 if (a->diag) { 2210 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2211 c->diag = a->diag; 2212 c->free_diag = PETSC_FALSE; 2213 } else { 2214 ierr = PetscMalloc(mbs*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2215 ierr = PetscLogObjectMemory(C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2216 for (i=0; i<mbs; i++) { 2217 c->diag[i] = a->diag[i]; 2218 } 2219 c->free_diag = PETSC_TRUE; 2220 } 2221 } else c->diag = 0; 2222 c->nz = a->nz; 2223 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2224 c->solve_work = 0; 2225 c->mult_work = 0; 2226 c->free_a = PETSC_TRUE; 2227 *B = C; 2228 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2229 PetscFunctionReturn(0); 2230 } 2231 2232 #undef __FUNCT__ 2233 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2234 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2235 { 2236 Mat_SeqSBAIJ *a; 2237 PetscErrorCode ierr; 2238 int fd; 2239 PetscMPIInt size; 2240 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2241 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2242 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2243 PetscInt *masked,nmask,tmp,bs2,ishift; 2244 PetscScalar *aa; 2245 MPI_Comm comm = ((PetscObject)viewer)->comm; 2246 2247 PetscFunctionBegin; 2248 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2249 bs2 = bs*bs; 2250 2251 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2252 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2253 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2254 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2255 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2256 M = header[1]; N = header[2]; nz = header[3]; 2257 2258 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2259 2260 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2261 2262 /* 2263 This code adds extra rows to make sure the number of rows is 2264 divisible by the blocksize 2265 */ 2266 mbs = M/bs; 2267 extra_rows = bs - M + bs*(mbs); 2268 if (extra_rows == bs) extra_rows = 0; 2269 else mbs++; 2270 if (extra_rows) { 2271 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2272 } 2273 2274 /* Set global sizes if not already set */ 2275 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2276 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2277 } else { /* Check if the matrix global sizes are correct */ 2278 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2279 if (M != rows || N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols); 2280 } 2281 2282 /* read in row lengths */ 2283 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2284 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2285 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2286 2287 /* read in column indices */ 2288 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2289 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2290 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2291 2292 /* loop over row lengths determining block row lengths */ 2293 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2294 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2295 ierr = PetscMalloc2(mbs,PetscInt,&mask,mbs,PetscInt,&masked);CHKERRQ(ierr); 2296 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2297 rowcount = 0; 2298 nzcount = 0; 2299 for (i=0; i<mbs; i++) { 2300 nmask = 0; 2301 for (j=0; j<bs; j++) { 2302 kmax = rowlengths[rowcount]; 2303 for (k=0; k<kmax; k++) { 2304 tmp = jj[nzcount++]/bs; /* block col. index */ 2305 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2306 } 2307 rowcount++; 2308 } 2309 s_browlengths[i] += nmask; 2310 2311 /* zero out the mask elements we set */ 2312 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2313 } 2314 2315 /* Do preallocation */ 2316 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2317 a = (Mat_SeqSBAIJ*)newmat->data; 2318 2319 /* set matrix "i" values */ 2320 a->i[0] = 0; 2321 for (i=1; i<= mbs; i++) { 2322 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2323 a->ilen[i-1] = s_browlengths[i-1]; 2324 } 2325 a->nz = a->i[mbs]; 2326 2327 /* read in nonzero values */ 2328 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2329 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2330 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2331 2332 /* set "a" and "j" values into matrix */ 2333 nzcount = 0; jcount = 0; 2334 for (i=0; i<mbs; i++) { 2335 nzcountb = nzcount; 2336 nmask = 0; 2337 for (j=0; j<bs; j++) { 2338 kmax = rowlengths[i*bs+j]; 2339 for (k=0; k<kmax; k++) { 2340 tmp = jj[nzcount++]/bs; /* block col. index */ 2341 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2342 } 2343 } 2344 /* sort the masked values */ 2345 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2346 2347 /* set "j" values into matrix */ 2348 maskcount = 1; 2349 for (j=0; j<nmask; j++) { 2350 a->j[jcount++] = masked[j]; 2351 mask[masked[j]] = maskcount++; 2352 } 2353 2354 /* set "a" values into matrix */ 2355 ishift = bs2*a->i[i]; 2356 for (j=0; j<bs; j++) { 2357 kmax = rowlengths[i*bs+j]; 2358 for (k=0; k<kmax; k++) { 2359 tmp = jj[nzcountb]/bs ; /* block col. index */ 2360 if (tmp >= i){ 2361 block = mask[tmp] - 1; 2362 point = jj[nzcountb] - bs*tmp; 2363 idx = ishift + bs2*block + j + bs*point; 2364 a->a[idx] = aa[nzcountb]; 2365 } 2366 nzcountb++; 2367 } 2368 } 2369 /* zero out the mask elements we set */ 2370 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2371 } 2372 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2373 2374 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2375 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2376 ierr = PetscFree(aa);CHKERRQ(ierr); 2377 ierr = PetscFree(jj);CHKERRQ(ierr); 2378 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2379 2380 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2381 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2382 ierr = MatView_Private(newmat);CHKERRQ(ierr); 2383 PetscFunctionReturn(0); 2384 } 2385 2386 #undef __FUNCT__ 2387 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2388 /*@ 2389 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2390 (upper triangular entries in CSR format) provided by the user. 2391 2392 Collective on MPI_Comm 2393 2394 Input Parameters: 2395 + comm - must be an MPI communicator of size 1 2396 . bs - size of block 2397 . m - number of rows 2398 . n - number of columns 2399 . i - row indices 2400 . j - column indices 2401 - a - matrix values 2402 2403 Output Parameter: 2404 . mat - the matrix 2405 2406 Level: advanced 2407 2408 Notes: 2409 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2410 once the matrix is destroyed 2411 2412 You cannot set new nonzero locations into this matrix, that will generate an error. 2413 2414 The i and j indices are 0 based 2415 2416 When block size is greater than 1 the matrix values must be stored using the SBAIJ storage format (see the SBAIJ code to determine this). For block size of 1 2417 it is the regular CSR format excluding the lower triangular elements. 2418 2419 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2420 2421 @*/ 2422 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2423 { 2424 PetscErrorCode ierr; 2425 PetscInt ii; 2426 Mat_SeqSBAIJ *sbaij; 2427 2428 PetscFunctionBegin; 2429 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2430 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2431 2432 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2433 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2434 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2435 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2436 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2437 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2438 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2439 2440 sbaij->i = i; 2441 sbaij->j = j; 2442 sbaij->a = a; 2443 sbaij->singlemalloc = PETSC_FALSE; 2444 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2445 sbaij->free_a = PETSC_FALSE; 2446 sbaij->free_ij = PETSC_FALSE; 2447 2448 for (ii=0; ii<m; ii++) { 2449 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2450 #if defined(PETSC_USE_DEBUG) 2451 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2452 #endif 2453 } 2454 #if defined(PETSC_USE_DEBUG) 2455 for (ii=0; ii<sbaij->i[m]; ii++) { 2456 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2457 if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2458 } 2459 #endif 2460 2461 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2462 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2463 PetscFunctionReturn(0); 2464 } 2465 2466 2467 2468 2469 2470