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]],PetscRealPart(1.0/a->a[diag[i]]),PetscImaginaryPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 414 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) { 415 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",a->j[diag[i]],PetscRealPart(1.0/a->a[diag[i]]),-PetscImaginaryPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 416 } else { 417 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[diag[i]],PetscRealPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 418 } 419 #else 420 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[diag[i]],1.0/a->a[diag[i]]);CHKERRQ(ierr); 421 #endif 422 /* off-diagonal entries */ 423 for (k=a->i[i]; k<a->i[i+1]-1; k++) { 424 #if defined(PETSC_USE_COMPLEX) 425 if (PetscImaginaryPart(a->a[k]) > 0.0) { 426 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k],PetscRealPart(a->a[k]),PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 427 } else if (PetscImaginaryPart(a->a[k]) < 0.0) { 428 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k],PetscRealPart(a->a[k]),-PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 429 } else { 430 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k],PetscRealPart(a->a[k]));CHKERRQ(ierr); 431 } 432 #else 433 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[k],a->a[k]);CHKERRQ(ierr); 434 #endif 435 } 436 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 437 } 438 439 } else { /* for non-factored matrix */ 440 for (i=0; i<a->mbs; i++) { /* for row block i */ 441 for (j=0; j<bs; j++) { /* for row bs*i + j */ 442 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 443 for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */ 444 for (l=0; l<bs; l++) { /* for column */ 445 #if defined(PETSC_USE_COMPLEX) 446 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 447 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 448 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 449 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 450 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 451 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 452 } else { 453 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 454 } 455 #else 456 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 457 #endif 458 } 459 } 460 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 461 } 462 } 463 } 464 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 465 } 466 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 467 PetscFunctionReturn(0); 468 } 469 470 #undef __FUNCT__ 471 #define __FUNCT__ "MatView_SeqSBAIJ_Draw_Zoom" 472 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 473 { 474 Mat A = (Mat) Aa; 475 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 476 PetscErrorCode ierr; 477 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 478 PetscMPIInt rank; 479 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 480 MatScalar *aa; 481 MPI_Comm comm; 482 PetscViewer viewer; 483 484 PetscFunctionBegin; 485 /* 486 This is nasty. If this is called from an originally parallel matrix 487 then all processes call this,but only the first has the matrix so the 488 rest should return immediately. 489 */ 490 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 491 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 492 if (rank) PetscFunctionReturn(0); 493 494 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 495 496 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 497 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 498 499 /* loop over matrix elements drawing boxes */ 500 color = PETSC_DRAW_BLUE; 501 for (i=0,row=0; i<mbs; i++,row+=bs) { 502 for (j=a->i[i]; j<a->i[i+1]; j++) { 503 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 504 x_l = a->j[j]*bs; x_r = x_l + 1.0; 505 aa = a->a + j*bs2; 506 for (k=0; k<bs; k++) { 507 for (l=0; l<bs; l++) { 508 if (PetscRealPart(*aa++) >= 0.) continue; 509 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 510 } 511 } 512 } 513 } 514 color = PETSC_DRAW_CYAN; 515 for (i=0,row=0; i<mbs; i++,row+=bs) { 516 for (j=a->i[i]; j<a->i[i+1]; j++) { 517 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 518 x_l = a->j[j]*bs; x_r = x_l + 1.0; 519 aa = a->a + j*bs2; 520 for (k=0; k<bs; k++) { 521 for (l=0; l<bs; l++) { 522 if (PetscRealPart(*aa++) != 0.) continue; 523 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 524 } 525 } 526 } 527 } 528 529 color = PETSC_DRAW_RED; 530 for (i=0,row=0; i<mbs; i++,row+=bs) { 531 for (j=a->i[i]; j<a->i[i+1]; j++) { 532 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 533 x_l = a->j[j]*bs; x_r = x_l + 1.0; 534 aa = a->a + j*bs2; 535 for (k=0; k<bs; k++) { 536 for (l=0; l<bs; l++) { 537 if (PetscRealPart(*aa++) <= 0.) continue; 538 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 539 } 540 } 541 } 542 } 543 PetscFunctionReturn(0); 544 } 545 546 #undef __FUNCT__ 547 #define __FUNCT__ "MatView_SeqSBAIJ_Draw" 548 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 549 { 550 PetscErrorCode ierr; 551 PetscReal xl,yl,xr,yr,w,h; 552 PetscDraw draw; 553 PetscBool isnull; 554 555 PetscFunctionBegin; 556 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 557 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 558 559 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 560 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 561 xr += w; yr += h; xl = -w; yl = -h; 562 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 563 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 564 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 565 PetscFunctionReturn(0); 566 } 567 568 #undef __FUNCT__ 569 #define __FUNCT__ "MatView_SeqSBAIJ" 570 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 571 { 572 PetscErrorCode ierr; 573 PetscBool iascii,isdraw; 574 FILE *file = 0; 575 576 PetscFunctionBegin; 577 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 578 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 579 if (iascii){ 580 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 581 } else if (isdraw) { 582 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 583 } else { 584 Mat B; 585 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 586 ierr = MatView(B,viewer);CHKERRQ(ierr); 587 ierr = MatDestroy(&B);CHKERRQ(ierr); 588 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 589 if (file) { 590 fprintf(file,"-matload_block_size %d\n",(int)A->rmap->bs); 591 } 592 } 593 PetscFunctionReturn(0); 594 } 595 596 597 #undef __FUNCT__ 598 #define __FUNCT__ "MatGetValues_SeqSBAIJ" 599 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 600 { 601 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 602 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 603 PetscInt *ai = a->i,*ailen = a->ilen; 604 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 605 MatScalar *ap,*aa = a->a; 606 607 PetscFunctionBegin; 608 for (k=0; k<m; k++) { /* loop over rows */ 609 row = im[k]; brow = row/bs; 610 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 611 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); 612 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 613 nrow = ailen[brow]; 614 for (l=0; l<n; l++) { /* loop over columns */ 615 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 616 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); 617 col = in[l] ; 618 bcol = col/bs; 619 cidx = col%bs; 620 ridx = row%bs; 621 high = nrow; 622 low = 0; /* assume unsorted */ 623 while (high-low > 5) { 624 t = (low+high)/2; 625 if (rp[t] > bcol) high = t; 626 else low = t; 627 } 628 for (i=low; i<high; i++) { 629 if (rp[i] > bcol) break; 630 if (rp[i] == bcol) { 631 *v++ = ap[bs2*i+bs*cidx+ridx]; 632 goto finished; 633 } 634 } 635 *v++ = 0.0; 636 finished:; 637 } 638 } 639 PetscFunctionReturn(0); 640 } 641 642 643 #undef __FUNCT__ 644 #define __FUNCT__ "MatSetValuesBlocked_SeqSBAIJ" 645 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 646 { 647 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 648 PetscErrorCode ierr; 649 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 650 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen; 651 PetscInt *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 652 PetscBool roworiented=a->roworiented; 653 const PetscScalar *value = v; 654 MatScalar *ap,*aa = a->a,*bap; 655 656 PetscFunctionBegin; 657 if (roworiented) { 658 stepval = (n-1)*bs; 659 } else { 660 stepval = (m-1)*bs; 661 } 662 for (k=0; k<m; k++) { /* loop over added rows */ 663 row = im[k]; 664 if (row < 0) continue; 665 #if defined(PETSC_USE_DEBUG) 666 if (row >= a->mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1); 667 #endif 668 rp = aj + ai[row]; 669 ap = aa + bs2*ai[row]; 670 rmax = imax[row]; 671 nrow = ailen[row]; 672 low = 0; 673 high = nrow; 674 for (l=0; l<n; l++) { /* loop over added columns */ 675 if (in[l] < 0) continue; 676 col = in[l]; 677 #if defined(PETSC_USE_DEBUG) 678 if (col >= a->nbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1); 679 #endif 680 if (col < row) { 681 if (a->ignore_ltriangular) { 682 continue; /* ignore lower triangular block */ 683 } else { 684 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)"); 685 } 686 } 687 if (roworiented) { 688 value = v + k*(stepval+bs)*bs + l*bs; 689 } else { 690 value = v + l*(stepval+bs)*bs + k*bs; 691 } 692 if (col <= lastcol) low = 0; else high = nrow; 693 lastcol = col; 694 while (high-low > 7) { 695 t = (low+high)/2; 696 if (rp[t] > col) high = t; 697 else low = t; 698 } 699 for (i=low; i<high; i++) { 700 if (rp[i] > col) break; 701 if (rp[i] == col) { 702 bap = ap + bs2*i; 703 if (roworiented) { 704 if (is == ADD_VALUES) { 705 for (ii=0; ii<bs; ii++,value+=stepval) { 706 for (jj=ii; jj<bs2; jj+=bs) { 707 bap[jj] += *value++; 708 } 709 } 710 } else { 711 for (ii=0; ii<bs; ii++,value+=stepval) { 712 for (jj=ii; jj<bs2; jj+=bs) { 713 bap[jj] = *value++; 714 } 715 } 716 } 717 } else { 718 if (is == ADD_VALUES) { 719 for (ii=0; ii<bs; ii++,value+=stepval) { 720 for (jj=0; jj<bs; jj++) { 721 *bap++ += *value++; 722 } 723 } 724 } else { 725 for (ii=0; ii<bs; ii++,value+=stepval) { 726 for (jj=0; jj<bs; jj++) { 727 *bap++ = *value++; 728 } 729 } 730 } 731 } 732 goto noinsert2; 733 } 734 } 735 if (nonew == 1) goto noinsert2; 736 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 737 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 738 N = nrow++ - 1; high++; 739 /* shift up all the later entries in this row */ 740 for (ii=N; ii>=i; ii--) { 741 rp[ii+1] = rp[ii]; 742 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 743 } 744 if (N >= i) { 745 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 746 } 747 rp[i] = col; 748 bap = ap + bs2*i; 749 if (roworiented) { 750 for (ii=0; ii<bs; ii++,value+=stepval) { 751 for (jj=ii; jj<bs2; jj+=bs) { 752 bap[jj] = *value++; 753 } 754 } 755 } else { 756 for (ii=0; ii<bs; ii++,value+=stepval) { 757 for (jj=0; jj<bs; jj++) { 758 *bap++ = *value++; 759 } 760 } 761 } 762 noinsert2:; 763 low = i; 764 } 765 ailen[row] = nrow; 766 } 767 PetscFunctionReturn(0); 768 } 769 770 /* 771 This is not yet used 772 */ 773 #undef __FUNCT__ 774 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode" 775 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 776 { 777 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 778 PetscErrorCode ierr; 779 const PetscInt *ai = a->i, *aj = a->j,*cols; 780 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 781 PetscBool flag; 782 783 PetscFunctionBegin; 784 ierr = PetscMalloc(m*sizeof(PetscInt),&ns);CHKERRQ(ierr); 785 while (i < m){ 786 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 787 /* Limits the number of elements in a node to 'a->inode.limit' */ 788 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 789 nzy = ai[j+1] - ai[j]; 790 if (nzy != (nzx - j + i)) break; 791 ierr = PetscMemcmp(aj + ai[i] + j - i,aj + ai[j],nzy*sizeof(PetscInt),&flag);CHKERRQ(ierr); 792 if (!flag) break; 793 } 794 ns[node_count++] = blk_size; 795 i = j; 796 } 797 if (!a->inode.size && m && node_count > .9*m) { 798 ierr = PetscFree(ns);CHKERRQ(ierr); 799 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 800 } else { 801 a->inode.node_count = node_count; 802 ierr = PetscMalloc(node_count*sizeof(PetscInt),&a->inode.size);CHKERRQ(ierr); 803 ierr = PetscLogObjectMemory(A,node_count*sizeof(PetscInt));CHKERRQ(ierr); 804 ierr = PetscMemcpy(a->inode.size,ns,node_count*sizeof(PetscInt)); 805 ierr = PetscFree(ns);CHKERRQ(ierr); 806 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 807 808 /* count collections of adjacent columns in each inode */ 809 row = 0; 810 cnt = 0; 811 for (i=0; i<node_count; i++) { 812 cols = aj + ai[row] + a->inode.size[i]; 813 nz = ai[row+1] - ai[row] - a->inode.size[i]; 814 for (j=1; j<nz; j++) { 815 if (cols[j] != cols[j-1]+1) { 816 cnt++; 817 } 818 } 819 cnt++; 820 row += a->inode.size[i]; 821 } 822 ierr = PetscMalloc(2*cnt*sizeof(PetscInt),&counts);CHKERRQ(ierr); 823 cnt = 0; 824 row = 0; 825 for (i=0; i<node_count; i++) { 826 cols = aj + ai[row] + a->inode.size[i]; 827 CHKMEMQ; 828 counts[2*cnt] = cols[0]; 829 CHKMEMQ; 830 nz = ai[row+1] - ai[row] - a->inode.size[i]; 831 cnt2 = 1; 832 for (j=1; j<nz; j++) { 833 if (cols[j] != cols[j-1]+1) { 834 CHKMEMQ; 835 counts[2*(cnt++)+1] = cnt2; 836 counts[2*cnt] = cols[j]; 837 CHKMEMQ; 838 cnt2 = 1; 839 } else cnt2++; 840 } 841 CHKMEMQ; 842 counts[2*(cnt++)+1] = cnt2; 843 CHKMEMQ; 844 row += a->inode.size[i]; 845 } 846 ierr = PetscIntView(2*cnt,counts,0); 847 } 848 PetscFunctionReturn(0); 849 } 850 851 #undef __FUNCT__ 852 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 853 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 854 { 855 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 856 PetscErrorCode ierr; 857 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 858 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 859 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 860 MatScalar *aa = a->a,*ap; 861 862 PetscFunctionBegin; 863 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 864 865 if (m) rmax = ailen[0]; 866 for (i=1; i<mbs; i++) { 867 /* move each row back by the amount of empty slots (fshift) before it*/ 868 fshift += imax[i-1] - ailen[i-1]; 869 rmax = PetscMax(rmax,ailen[i]); 870 if (fshift) { 871 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 872 N = ailen[i]; 873 for (j=0; j<N; j++) { 874 ip[j-fshift] = ip[j]; 875 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 876 } 877 } 878 ai[i] = ai[i-1] + ailen[i-1]; 879 } 880 if (mbs) { 881 fshift += imax[mbs-1] - ailen[mbs-1]; 882 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 883 } 884 /* reset ilen and imax for each row */ 885 for (i=0; i<mbs; i++) { 886 ailen[i] = imax[i] = ai[i+1] - ai[i]; 887 } 888 a->nz = ai[mbs]; 889 890 /* diagonals may have moved, reset it */ 891 if (a->diag) { 892 ierr = PetscMemcpy(a->diag,ai,mbs*sizeof(PetscInt));CHKERRQ(ierr); 893 } 894 if (fshift && a->nounused == -1) { 895 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); 896 } 897 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); 898 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 899 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 900 A->info.mallocs += a->reallocs; 901 a->reallocs = 0; 902 A->info.nz_unneeded = (PetscReal)fshift*bs2; 903 a->idiagvalid = PETSC_FALSE; 904 905 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 906 if (!a->jshort) { 907 ierr = PetscMalloc(a->i[A->rmap->n]*sizeof(unsigned short),&a->jshort);CHKERRQ(ierr); 908 ierr = PetscLogObjectMemory(A,a->i[A->rmap->n]*sizeof(unsigned short));CHKERRQ(ierr); 909 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 910 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 911 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 912 a->free_jshort = PETSC_TRUE; 913 } 914 } 915 PetscFunctionReturn(0); 916 } 917 918 /* 919 This function returns an array of flags which indicate the locations of contiguous 920 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 921 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 922 Assume: sizes should be long enough to hold all the values. 923 */ 924 #undef __FUNCT__ 925 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 926 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 927 { 928 PetscInt i,j,k,row; 929 PetscBool flg; 930 931 PetscFunctionBegin; 932 for (i=0,j=0; i<n; j++) { 933 row = idx[i]; 934 if (row%bs!=0) { /* Not the begining of a block */ 935 sizes[j] = 1; 936 i++; 937 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 938 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 939 i++; 940 } else { /* Begining of the block, so check if the complete block exists */ 941 flg = PETSC_TRUE; 942 for (k=1; k<bs; k++) { 943 if (row+k != idx[i+k]) { /* break in the block */ 944 flg = PETSC_FALSE; 945 break; 946 } 947 } 948 if (flg) { /* No break in the bs */ 949 sizes[j] = bs; 950 i+= bs; 951 } else { 952 sizes[j] = 1; 953 i++; 954 } 955 } 956 } 957 *bs_max = j; 958 PetscFunctionReturn(0); 959 } 960 961 962 /* Only add/insert a(i,j) with i<=j (blocks). 963 Any a(i,j) with i>j input by user is ingored. 964 */ 965 966 #undef __FUNCT__ 967 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 968 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 969 { 970 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 971 PetscErrorCode ierr; 972 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 973 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 974 PetscInt *aj=a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 975 PetscInt ridx,cidx,bs2=a->bs2; 976 MatScalar *ap,value,*aa=a->a,*bap; 977 978 PetscFunctionBegin; 979 if (v) PetscValidScalarPointer(v,6); 980 for (k=0; k<m; k++) { /* loop over added rows */ 981 row = im[k]; /* row number */ 982 brow = row/bs; /* block row number */ 983 if (row < 0) continue; 984 #if defined(PETSC_USE_DEBUG) 985 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); 986 #endif 987 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 988 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 989 rmax = imax[brow]; /* maximum space allocated for this row */ 990 nrow = ailen[brow]; /* actual length of this row */ 991 low = 0; 992 993 for (l=0; l<n; l++) { /* loop over added columns */ 994 if (in[l] < 0) continue; 995 #if defined(PETSC_USE_DEBUG) 996 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); 997 #endif 998 col = in[l]; 999 bcol = col/bs; /* block col number */ 1000 1001 if (brow > bcol) { 1002 if (a->ignore_ltriangular){ 1003 continue; /* ignore lower triangular values */ 1004 } else { 1005 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)"); 1006 } 1007 } 1008 1009 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 1010 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 1011 /* element value a(k,l) */ 1012 if (roworiented) { 1013 value = v[l + k*n]; 1014 } else { 1015 value = v[k + l*m]; 1016 } 1017 1018 /* move pointer bap to a(k,l) quickly and add/insert value */ 1019 if (col <= lastcol) low = 0; high = nrow; 1020 lastcol = col; 1021 while (high-low > 7) { 1022 t = (low+high)/2; 1023 if (rp[t] > bcol) high = t; 1024 else low = t; 1025 } 1026 for (i=low; i<high; i++) { 1027 if (rp[i] > bcol) break; 1028 if (rp[i] == bcol) { 1029 bap = ap + bs2*i + bs*cidx + ridx; 1030 if (is == ADD_VALUES) *bap += value; 1031 else *bap = value; 1032 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1033 if (brow == bcol && ridx < cidx){ 1034 bap = ap + bs2*i + bs*ridx + cidx; 1035 if (is == ADD_VALUES) *bap += value; 1036 else *bap = value; 1037 } 1038 goto noinsert1; 1039 } 1040 } 1041 1042 if (nonew == 1) goto noinsert1; 1043 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 1044 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 1045 1046 N = nrow++ - 1; high++; 1047 /* shift up all the later entries in this row */ 1048 for (ii=N; ii>=i; ii--) { 1049 rp[ii+1] = rp[ii]; 1050 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 1051 } 1052 if (N>=i) { 1053 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1054 } 1055 rp[i] = bcol; 1056 ap[bs2*i + bs*cidx + ridx] = value; 1057 noinsert1:; 1058 low = i; 1059 } 1060 } /* end of loop over added columns */ 1061 ailen[brow] = nrow; 1062 } /* end of loop over added rows */ 1063 PetscFunctionReturn(0); 1064 } 1065 1066 #undef __FUNCT__ 1067 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 1068 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1069 { 1070 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1071 Mat outA; 1072 PetscErrorCode ierr; 1073 PetscBool row_identity; 1074 1075 PetscFunctionBegin; 1076 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1077 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1078 if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1079 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()! */ 1080 1081 outA = inA; 1082 inA->factortype = MAT_FACTOR_ICC; 1083 1084 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1085 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr); 1086 1087 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1088 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1089 a->row = row; 1090 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1091 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1092 a->col = row; 1093 1094 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1095 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1096 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1097 1098 if (!a->solve_work) { 1099 ierr = PetscMalloc((inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1100 ierr = PetscLogObjectMemory(inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1101 } 1102 1103 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 EXTERN_C_BEGIN 1108 #undef __FUNCT__ 1109 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1110 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1111 { 1112 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1113 PetscInt i,nz,n; 1114 PetscErrorCode ierr; 1115 1116 PetscFunctionBegin; 1117 nz = baij->maxnz; 1118 n = mat->cmap->n; 1119 for (i=0; i<nz; i++) { 1120 baij->j[i] = indices[i]; 1121 } 1122 baij->nz = nz; 1123 for (i=0; i<n; i++) { 1124 baij->ilen[i] = baij->imax[i]; 1125 } 1126 ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1127 PetscFunctionReturn(0); 1128 } 1129 EXTERN_C_END 1130 1131 #undef __FUNCT__ 1132 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1133 /*@ 1134 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1135 in the matrix. 1136 1137 Input Parameters: 1138 + mat - the SeqSBAIJ matrix 1139 - indices - the column indices 1140 1141 Level: advanced 1142 1143 Notes: 1144 This can be called if you have precomputed the nonzero structure of the 1145 matrix and want to provide it to the matrix object to improve the performance 1146 of the MatSetValues() operation. 1147 1148 You MUST have set the correct numbers of nonzeros per row in the call to 1149 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1150 1151 MUST be called before any calls to MatSetValues() 1152 1153 .seealso: MatCreateSeqSBAIJ 1154 @*/ 1155 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1156 { 1157 PetscErrorCode ierr; 1158 1159 PetscFunctionBegin; 1160 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1161 PetscValidPointer(indices,2); 1162 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1168 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1169 { 1170 PetscErrorCode ierr; 1171 1172 PetscFunctionBegin; 1173 /* If the two matrices have the same copy implementation, use fast copy. */ 1174 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1175 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1176 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1177 1178 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"); 1179 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1180 } else { 1181 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1182 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1183 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1184 } 1185 PetscFunctionReturn(0); 1186 } 1187 1188 #undef __FUNCT__ 1189 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1190 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1191 { 1192 PetscErrorCode ierr; 1193 1194 PetscFunctionBegin; 1195 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,-PetscMax(A->rmap->bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1196 PetscFunctionReturn(0); 1197 } 1198 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1201 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1202 { 1203 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1204 PetscFunctionBegin; 1205 *array = a->a; 1206 PetscFunctionReturn(0); 1207 } 1208 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1211 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1212 { 1213 PetscFunctionBegin; 1214 PetscFunctionReturn(0); 1215 } 1216 1217 #undef __FUNCT__ 1218 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1219 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1220 { 1221 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1222 PetscErrorCode ierr; 1223 PetscInt i,bs=Y->rmap->bs,bs2=bs*bs,j; 1224 PetscBLASInt one = 1; 1225 1226 PetscFunctionBegin; 1227 if (str == SAME_NONZERO_PATTERN) { 1228 PetscScalar alpha = a; 1229 PetscBLASInt bnz = PetscBLASIntCast(x->nz*bs2); 1230 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1231 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1232 if (y->xtoy && y->XtoY != X) { 1233 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1234 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 1235 } 1236 if (!y->xtoy) { /* get xtoy */ 1237 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1238 y->XtoY = X; 1239 } 1240 for (i=0; i<x->nz; i++) { 1241 j = 0; 1242 while (j < bs2){ 1243 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1244 j++; 1245 } 1246 } 1247 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); 1248 } else { 1249 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1250 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1251 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1252 } 1253 PetscFunctionReturn(0); 1254 } 1255 1256 #undef __FUNCT__ 1257 #define __FUNCT__ "MatSetBlockSize_SeqSBAIJ" 1258 PetscErrorCode MatSetBlockSize_SeqSBAIJ(Mat A,PetscInt bs) 1259 { 1260 PetscInt rbs,cbs; 1261 PetscErrorCode ierr; 1262 1263 PetscFunctionBegin; 1264 ierr = PetscLayoutGetBlockSize(A->rmap,&rbs);CHKERRQ(ierr); 1265 ierr = PetscLayoutGetBlockSize(A->cmap,&cbs);CHKERRQ(ierr); 1266 if (rbs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Attempt to set block size %d with SBAIJ %d",bs,rbs); 1267 if (cbs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Attempt to set block size %d with SBAIJ %d",bs,cbs); 1268 PetscFunctionReturn(0); 1269 } 1270 1271 #undef __FUNCT__ 1272 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1273 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1274 { 1275 PetscFunctionBegin; 1276 *flg = PETSC_TRUE; 1277 PetscFunctionReturn(0); 1278 } 1279 1280 #undef __FUNCT__ 1281 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1282 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1283 { 1284 PetscFunctionBegin; 1285 *flg = PETSC_TRUE; 1286 PetscFunctionReturn(0); 1287 } 1288 1289 #undef __FUNCT__ 1290 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1291 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1292 { 1293 PetscFunctionBegin; 1294 *flg = PETSC_FALSE; 1295 PetscFunctionReturn(0); 1296 } 1297 1298 #undef __FUNCT__ 1299 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1300 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1301 { 1302 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1303 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1304 MatScalar *aa = a->a; 1305 1306 PetscFunctionBegin; 1307 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1308 PetscFunctionReturn(0); 1309 } 1310 1311 #undef __FUNCT__ 1312 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1313 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1314 { 1315 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1316 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1317 MatScalar *aa = a->a; 1318 1319 PetscFunctionBegin; 1320 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1321 PetscFunctionReturn(0); 1322 } 1323 1324 #undef __FUNCT__ 1325 #define __FUNCT__ "MatZeroRowsColumns_SeqSBAIJ" 1326 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1327 { 1328 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1329 PetscErrorCode ierr; 1330 PetscInt i,j,k,count; 1331 PetscInt bs=A->rmap->bs,bs2=baij->bs2,row,col; 1332 PetscScalar zero = 0.0; 1333 MatScalar *aa; 1334 const PetscScalar *xx; 1335 PetscScalar *bb; 1336 PetscBool *zeroed,vecs = PETSC_FALSE; 1337 1338 PetscFunctionBegin; 1339 /* fix right hand side if needed */ 1340 if (x && b) { 1341 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1342 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1343 vecs = PETSC_TRUE; 1344 } 1345 A->same_nonzero = PETSC_TRUE; 1346 1347 /* zero the columns */ 1348 ierr = PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);CHKERRQ(ierr); 1349 ierr = PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));CHKERRQ(ierr); 1350 for (i=0; i<is_n; i++) { 1351 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]); 1352 zeroed[is_idx[i]] = PETSC_TRUE; 1353 } 1354 if (vecs) { 1355 for (i=0; i<A->rmap->N; i++) { 1356 row = i/bs; 1357 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1358 for (k=0; k<bs; k++) { 1359 col = bs*baij->j[j] + k; 1360 if (col <= i) continue; 1361 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1362 if (!zeroed[i] && zeroed[col]) { 1363 bb[i] -= aa[0]*xx[col]; 1364 } 1365 if (zeroed[i] && !zeroed[col]) { 1366 bb[col] -= aa[0]*xx[i]; 1367 } 1368 } 1369 } 1370 } 1371 for (i=0; i<is_n; i++) { 1372 bb[is_idx[i]] = diag*xx[is_idx[i]]; 1373 } 1374 } 1375 1376 for (i=0; i<A->rmap->N; i++) { 1377 if (!zeroed[i]) { 1378 row = i/bs; 1379 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1380 for (k=0; k<bs; k++) { 1381 col = bs*baij->j[j] + k; 1382 if (zeroed[col]) { 1383 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1384 aa[0] = 0.0; 1385 } 1386 } 1387 } 1388 } 1389 } 1390 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1391 if (vecs) { 1392 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1393 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1394 } 1395 1396 /* zero the rows */ 1397 for (i=0; i<is_n; i++) { 1398 row = is_idx[i]; 1399 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1400 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1401 for (k=0; k<count; k++) { 1402 aa[0] = zero; 1403 aa += bs; 1404 } 1405 if (diag != 0.0) { 1406 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1407 } 1408 } 1409 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1410 PetscFunctionReturn(0); 1411 } 1412 1413 /* -------------------------------------------------------------------*/ 1414 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1415 MatGetRow_SeqSBAIJ, 1416 MatRestoreRow_SeqSBAIJ, 1417 MatMult_SeqSBAIJ_N, 1418 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1419 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1420 MatMultAdd_SeqSBAIJ_N, 1421 0, 1422 0, 1423 0, 1424 /*10*/ 0, 1425 0, 1426 MatCholeskyFactor_SeqSBAIJ, 1427 MatSOR_SeqSBAIJ, 1428 MatTranspose_SeqSBAIJ, 1429 /*15*/ MatGetInfo_SeqSBAIJ, 1430 MatEqual_SeqSBAIJ, 1431 MatGetDiagonal_SeqSBAIJ, 1432 MatDiagonalScale_SeqSBAIJ, 1433 MatNorm_SeqSBAIJ, 1434 /*20*/ 0, 1435 MatAssemblyEnd_SeqSBAIJ, 1436 MatSetOption_SeqSBAIJ, 1437 MatZeroEntries_SeqSBAIJ, 1438 /*24*/ 0, 1439 0, 1440 0, 1441 0, 1442 0, 1443 /*29*/ MatSetUpPreallocation_SeqSBAIJ, 1444 0, 1445 0, 1446 MatGetArray_SeqSBAIJ, 1447 MatRestoreArray_SeqSBAIJ, 1448 /*34*/ MatDuplicate_SeqSBAIJ, 1449 0, 1450 0, 1451 0, 1452 MatICCFactor_SeqSBAIJ, 1453 /*39*/ MatAXPY_SeqSBAIJ, 1454 MatGetSubMatrices_SeqSBAIJ, 1455 MatIncreaseOverlap_SeqSBAIJ, 1456 MatGetValues_SeqSBAIJ, 1457 MatCopy_SeqSBAIJ, 1458 /*44*/ 0, 1459 MatScale_SeqSBAIJ, 1460 0, 1461 0, 1462 MatZeroRowsColumns_SeqSBAIJ, 1463 /*49*/ MatSetBlockSize_SeqSBAIJ, 1464 MatGetRowIJ_SeqSBAIJ, 1465 MatRestoreRowIJ_SeqSBAIJ, 1466 0, 1467 0, 1468 /*54*/ 0, 1469 0, 1470 0, 1471 0, 1472 MatSetValuesBlocked_SeqSBAIJ, 1473 /*59*/ MatGetSubMatrix_SeqSBAIJ, 1474 0, 1475 0, 1476 0, 1477 0, 1478 /*64*/ 0, 1479 0, 1480 0, 1481 0, 1482 0, 1483 /*69*/ MatGetRowMaxAbs_SeqSBAIJ, 1484 0, 1485 0, 1486 0, 1487 0, 1488 /*74*/ 0, 1489 0, 1490 0, 1491 0, 1492 0, 1493 /*79*/ 0, 1494 0, 1495 0, 1496 MatGetInertia_SeqSBAIJ, 1497 MatLoad_SeqSBAIJ, 1498 /*84*/ MatIsSymmetric_SeqSBAIJ, 1499 MatIsHermitian_SeqSBAIJ, 1500 MatIsStructurallySymmetric_SeqSBAIJ, 1501 0, 1502 0, 1503 /*89*/ 0, 1504 0, 1505 0, 1506 0, 1507 0, 1508 /*94*/ 0, 1509 0, 1510 0, 1511 0, 1512 0, 1513 /*99*/ 0, 1514 0, 1515 0, 1516 0, 1517 0, 1518 /*104*/0, 1519 MatRealPart_SeqSBAIJ, 1520 MatImaginaryPart_SeqSBAIJ, 1521 MatGetRowUpperTriangular_SeqSBAIJ, 1522 MatRestoreRowUpperTriangular_SeqSBAIJ, 1523 /*109*/0, 1524 0, 1525 0, 1526 0, 1527 MatMissingDiagonal_SeqSBAIJ, 1528 /*114*/0, 1529 0, 1530 0, 1531 0, 1532 0, 1533 /*119*/0, 1534 0, 1535 0, 1536 0 1537 }; 1538 1539 EXTERN_C_BEGIN 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1542 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1543 { 1544 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1545 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1546 PetscErrorCode ierr; 1547 1548 PetscFunctionBegin; 1549 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1550 1551 /* allocate space for values if not already there */ 1552 if (!aij->saved_values) { 1553 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1554 } 1555 1556 /* copy values over */ 1557 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1558 PetscFunctionReturn(0); 1559 } 1560 EXTERN_C_END 1561 1562 EXTERN_C_BEGIN 1563 #undef __FUNCT__ 1564 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1565 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1566 { 1567 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1568 PetscErrorCode ierr; 1569 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1570 1571 PetscFunctionBegin; 1572 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1573 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1574 1575 /* copy values over */ 1576 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1577 PetscFunctionReturn(0); 1578 } 1579 EXTERN_C_END 1580 1581 EXTERN_C_BEGIN 1582 #undef __FUNCT__ 1583 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1584 PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1585 { 1586 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1587 PetscErrorCode ierr; 1588 PetscInt i,mbs,bs2, newbs = PetscAbs(bs); 1589 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1590 1591 PetscFunctionBegin; 1592 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1593 B->preallocated = PETSC_TRUE; 1594 if (bs < 0) { 1595 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for MPISBAIJ matrix","Mat");CHKERRQ(ierr); 1596 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSeqSBAIJSetPreallocation",newbs,&newbs,PETSC_NULL);CHKERRQ(ierr); 1597 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1598 bs = PetscAbs(bs); 1599 } 1600 if (nnz && newbs != bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot change blocksize from command line if setting nnz"); 1601 bs = newbs; 1602 1603 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1604 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1605 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1606 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1607 1608 mbs = B->rmap->N/bs; 1609 bs2 = bs*bs; 1610 1611 if (mbs*bs != B->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1612 1613 if (nz == MAT_SKIP_ALLOCATION) { 1614 skipallocation = PETSC_TRUE; 1615 nz = 0; 1616 } 1617 1618 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1619 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1620 if (nnz) { 1621 for (i=0; i<mbs; i++) { 1622 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]); 1623 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); 1624 } 1625 } 1626 1627 B->ops->mult = MatMult_SeqSBAIJ_N; 1628 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1629 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1630 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1631 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1632 if (!flg) { 1633 switch (bs) { 1634 case 1: 1635 B->ops->mult = MatMult_SeqSBAIJ_1; 1636 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1637 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1638 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1639 break; 1640 case 2: 1641 B->ops->mult = MatMult_SeqSBAIJ_2; 1642 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1643 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1644 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1645 break; 1646 case 3: 1647 B->ops->mult = MatMult_SeqSBAIJ_3; 1648 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1649 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1650 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1651 break; 1652 case 4: 1653 B->ops->mult = MatMult_SeqSBAIJ_4; 1654 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1655 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1656 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1657 break; 1658 case 5: 1659 B->ops->mult = MatMult_SeqSBAIJ_5; 1660 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1661 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1662 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1663 break; 1664 case 6: 1665 B->ops->mult = MatMult_SeqSBAIJ_6; 1666 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1667 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1668 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1669 break; 1670 case 7: 1671 B->ops->mult = MatMult_SeqSBAIJ_7; 1672 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1673 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1674 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1675 break; 1676 } 1677 } 1678 1679 b->mbs = mbs; 1680 b->nbs = mbs; 1681 if (!skipallocation) { 1682 if (!b->imax) { 1683 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1684 b->free_imax_ilen = PETSC_TRUE; 1685 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1686 } 1687 if (!nnz) { 1688 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1689 else if (nz <= 0) nz = 1; 1690 for (i=0; i<mbs; i++) { 1691 b->imax[i] = nz; 1692 } 1693 nz = nz*mbs; /* total nz */ 1694 } else { 1695 nz = 0; 1696 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1697 } 1698 /* b->ilen will count nonzeros in each block row so far. */ 1699 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1700 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1701 1702 /* allocate the matrix space */ 1703 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1704 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1705 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1706 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1707 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1708 b->singlemalloc = PETSC_TRUE; 1709 1710 /* pointer to beginning of each row */ 1711 b->i[0] = 0; 1712 for (i=1; i<mbs+1; i++) { 1713 b->i[i] = b->i[i-1] + b->imax[i-1]; 1714 } 1715 b->free_a = PETSC_TRUE; 1716 b->free_ij = PETSC_TRUE; 1717 } else { 1718 b->free_a = PETSC_FALSE; 1719 b->free_ij = PETSC_FALSE; 1720 } 1721 1722 B->rmap->bs = bs; 1723 b->bs2 = bs2; 1724 b->nz = 0; 1725 b->maxnz = nz; 1726 1727 b->inew = 0; 1728 b->jnew = 0; 1729 b->anew = 0; 1730 b->a2anew = 0; 1731 b->permute = PETSC_FALSE; 1732 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 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 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2060 PetscValidType(B,1); 2061 PetscValidLogicalCollectiveInt(B,bs,2); 2062 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2063 PetscFunctionReturn(0); 2064 } 2065 2066 #undef __FUNCT__ 2067 #define __FUNCT__ "MatCreateSeqSBAIJ" 2068 /*@C 2069 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2070 compressed row) format. For good matrix assembly performance the 2071 user should preallocate the matrix storage by setting the parameter nz 2072 (or the array nnz). By setting these parameters accurately, performance 2073 during matrix assembly can be increased by more than a factor of 50. 2074 2075 Collective on MPI_Comm 2076 2077 Input Parameters: 2078 + comm - MPI communicator, set to PETSC_COMM_SELF 2079 . bs - size of block 2080 . m - number of rows, or number of columns 2081 . nz - number of block nonzeros per block row (same for all rows) 2082 - nnz - array containing the number of block nonzeros in the upper triangular plus 2083 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 2084 2085 Output Parameter: 2086 . A - the symmetric matrix 2087 2088 Options Database Keys: 2089 . -mat_no_unroll - uses code that does not unroll the loops in the 2090 block calculations (much slower) 2091 . -mat_block_size - size of the blocks to use 2092 2093 Level: intermediate 2094 2095 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2096 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2097 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2098 2099 Notes: 2100 The number of rows and columns must be divisible by blocksize. 2101 This matrix type does not support complex Hermitian operation. 2102 2103 Specify the preallocated storage with either nz or nnz (not both). 2104 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2105 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2106 2107 If the nnz parameter is given then the nz parameter is ignored 2108 2109 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 2110 @*/ 2111 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2112 { 2113 PetscErrorCode ierr; 2114 2115 PetscFunctionBegin; 2116 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2117 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2118 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2119 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2120 PetscFunctionReturn(0); 2121 } 2122 2123 #undef __FUNCT__ 2124 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2125 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2126 { 2127 Mat C; 2128 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2129 PetscErrorCode ierr; 2130 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2131 2132 PetscFunctionBegin; 2133 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2134 2135 *B = 0; 2136 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 2137 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2138 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2139 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2140 c = (Mat_SeqSBAIJ*)C->data; 2141 2142 C->preallocated = PETSC_TRUE; 2143 C->factortype = A->factortype; 2144 c->row = 0; 2145 c->icol = 0; 2146 c->saved_values = 0; 2147 c->keepnonzeropattern = a->keepnonzeropattern; 2148 C->assembled = PETSC_TRUE; 2149 2150 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2151 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2152 c->bs2 = a->bs2; 2153 c->mbs = a->mbs; 2154 c->nbs = a->nbs; 2155 2156 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2157 c->imax = a->imax; 2158 c->ilen = a->ilen; 2159 c->free_imax_ilen = PETSC_FALSE; 2160 } else { 2161 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 2162 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2163 for (i=0; i<mbs; i++) { 2164 c->imax[i] = a->imax[i]; 2165 c->ilen[i] = a->ilen[i]; 2166 } 2167 c->free_imax_ilen = PETSC_TRUE; 2168 } 2169 2170 /* allocate the matrix space */ 2171 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2172 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2173 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2174 c->singlemalloc = PETSC_FALSE; 2175 c->free_ij = PETSC_FALSE; 2176 c->parent = A; 2177 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2178 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2179 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2180 } else { 2181 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2182 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2183 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2184 c->singlemalloc = PETSC_TRUE; 2185 c->free_ij = PETSC_TRUE; 2186 } 2187 if (mbs > 0) { 2188 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2189 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2190 } 2191 if (cpvalues == MAT_COPY_VALUES) { 2192 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2193 } else { 2194 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2195 } 2196 if (a->jshort) { 2197 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2198 c->jshort = a->jshort; 2199 c->free_jshort = PETSC_FALSE; 2200 } else { 2201 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2202 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2203 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2204 c->free_jshort = PETSC_TRUE; 2205 } 2206 } 2207 } 2208 2209 c->roworiented = a->roworiented; 2210 c->nonew = a->nonew; 2211 2212 if (a->diag) { 2213 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2214 c->diag = a->diag; 2215 c->free_diag = PETSC_FALSE; 2216 } else { 2217 ierr = PetscMalloc(mbs*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2218 ierr = PetscLogObjectMemory(C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2219 for (i=0; i<mbs; i++) { 2220 c->diag[i] = a->diag[i]; 2221 } 2222 c->free_diag = PETSC_TRUE; 2223 } 2224 } else c->diag = 0; 2225 c->nz = a->nz; 2226 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2227 c->solve_work = 0; 2228 c->mult_work = 0; 2229 c->free_a = PETSC_TRUE; 2230 *B = C; 2231 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2232 PetscFunctionReturn(0); 2233 } 2234 2235 #undef __FUNCT__ 2236 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2237 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2238 { 2239 Mat_SeqSBAIJ *a; 2240 PetscErrorCode ierr; 2241 int fd; 2242 PetscMPIInt size; 2243 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2244 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2245 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2246 PetscInt *masked,nmask,tmp,bs2,ishift; 2247 PetscScalar *aa; 2248 MPI_Comm comm = ((PetscObject)viewer)->comm; 2249 2250 PetscFunctionBegin; 2251 ierr = PetscOptionsGetInt(((PetscObject)newmat)->prefix,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2252 bs2 = bs*bs; 2253 2254 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2255 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2256 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2257 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2258 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2259 M = header[1]; N = header[2]; nz = header[3]; 2260 2261 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2262 2263 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2264 2265 /* 2266 This code adds extra rows to make sure the number of rows is 2267 divisible by the blocksize 2268 */ 2269 mbs = M/bs; 2270 extra_rows = bs - M + bs*(mbs); 2271 if (extra_rows == bs) extra_rows = 0; 2272 else mbs++; 2273 if (extra_rows) { 2274 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2275 } 2276 2277 /* Set global sizes if not already set */ 2278 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2279 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2280 } else { /* Check if the matrix global sizes are correct */ 2281 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2282 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); 2283 } 2284 2285 /* read in row lengths */ 2286 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2287 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2288 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2289 2290 /* read in column indices */ 2291 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2292 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2293 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2294 2295 /* loop over row lengths determining block row lengths */ 2296 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2297 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2298 ierr = PetscMalloc2(mbs,PetscInt,&mask,mbs,PetscInt,&masked);CHKERRQ(ierr); 2299 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2300 rowcount = 0; 2301 nzcount = 0; 2302 for (i=0; i<mbs; i++) { 2303 nmask = 0; 2304 for (j=0; j<bs; j++) { 2305 kmax = rowlengths[rowcount]; 2306 for (k=0; k<kmax; k++) { 2307 tmp = jj[nzcount++]/bs; /* block col. index */ 2308 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2309 } 2310 rowcount++; 2311 } 2312 s_browlengths[i] += nmask; 2313 2314 /* zero out the mask elements we set */ 2315 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2316 } 2317 2318 /* Do preallocation */ 2319 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2320 a = (Mat_SeqSBAIJ*)newmat->data; 2321 2322 /* set matrix "i" values */ 2323 a->i[0] = 0; 2324 for (i=1; i<= mbs; i++) { 2325 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2326 a->ilen[i-1] = s_browlengths[i-1]; 2327 } 2328 a->nz = a->i[mbs]; 2329 2330 /* read in nonzero values */ 2331 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2332 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2333 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2334 2335 /* set "a" and "j" values into matrix */ 2336 nzcount = 0; jcount = 0; 2337 for (i=0; i<mbs; i++) { 2338 nzcountb = nzcount; 2339 nmask = 0; 2340 for (j=0; j<bs; j++) { 2341 kmax = rowlengths[i*bs+j]; 2342 for (k=0; k<kmax; k++) { 2343 tmp = jj[nzcount++]/bs; /* block col. index */ 2344 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2345 } 2346 } 2347 /* sort the masked values */ 2348 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2349 2350 /* set "j" values into matrix */ 2351 maskcount = 1; 2352 for (j=0; j<nmask; j++) { 2353 a->j[jcount++] = masked[j]; 2354 mask[masked[j]] = maskcount++; 2355 } 2356 2357 /* set "a" values into matrix */ 2358 ishift = bs2*a->i[i]; 2359 for (j=0; j<bs; j++) { 2360 kmax = rowlengths[i*bs+j]; 2361 for (k=0; k<kmax; k++) { 2362 tmp = jj[nzcountb]/bs ; /* block col. index */ 2363 if (tmp >= i){ 2364 block = mask[tmp] - 1; 2365 point = jj[nzcountb] - bs*tmp; 2366 idx = ishift + bs2*block + j + bs*point; 2367 a->a[idx] = aa[nzcountb]; 2368 } 2369 nzcountb++; 2370 } 2371 } 2372 /* zero out the mask elements we set */ 2373 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2374 } 2375 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2376 2377 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2378 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2379 ierr = PetscFree(aa);CHKERRQ(ierr); 2380 ierr = PetscFree(jj);CHKERRQ(ierr); 2381 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2382 2383 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2384 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2385 ierr = MatView_Private(newmat);CHKERRQ(ierr); 2386 PetscFunctionReturn(0); 2387 } 2388 2389 #undef __FUNCT__ 2390 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2391 /*@ 2392 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2393 (upper triangular entries in CSR format) provided by the user. 2394 2395 Collective on MPI_Comm 2396 2397 Input Parameters: 2398 + comm - must be an MPI communicator of size 1 2399 . bs - size of block 2400 . m - number of rows 2401 . n - number of columns 2402 . i - row indices 2403 . j - column indices 2404 - a - matrix values 2405 2406 Output Parameter: 2407 . mat - the matrix 2408 2409 Level: advanced 2410 2411 Notes: 2412 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2413 once the matrix is destroyed 2414 2415 You cannot set new nonzero locations into this matrix, that will generate an error. 2416 2417 The i and j indices are 0 based 2418 2419 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 2420 it is the regular CSR format excluding the lower triangular elements. 2421 2422 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2423 2424 @*/ 2425 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2426 { 2427 PetscErrorCode ierr; 2428 PetscInt ii; 2429 Mat_SeqSBAIJ *sbaij; 2430 2431 PetscFunctionBegin; 2432 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2433 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2434 2435 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2436 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2437 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2438 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2439 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2440 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2441 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2442 2443 sbaij->i = i; 2444 sbaij->j = j; 2445 sbaij->a = a; 2446 sbaij->singlemalloc = PETSC_FALSE; 2447 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2448 sbaij->free_a = PETSC_FALSE; 2449 sbaij->free_ij = PETSC_FALSE; 2450 2451 for (ii=0; ii<m; ii++) { 2452 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2453 #if defined(PETSC_USE_DEBUG) 2454 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]); 2455 #endif 2456 } 2457 #if defined(PETSC_USE_DEBUG) 2458 for (ii=0; ii<sbaij->i[m]; ii++) { 2459 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2460 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]); 2461 } 2462 #endif 2463 2464 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2465 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2466 PetscFunctionReturn(0); 2467 } 2468 2469 2470 2471 2472 2473