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