1 /* 2 Defines the basic matrix operations for the AIJ (compressed row) 3 matrix storage format. 4 */ 5 6 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 #include <petscbt.h> 9 #include <petsc/private/kernels/blocktranspose.h> 10 11 PetscErrorCode MatSeqAIJSetTypeFromOptions(Mat A) 12 { 13 PetscBool flg; 14 char type[256]; 15 16 PetscFunctionBegin; 17 PetscObjectOptionsBegin((PetscObject)A); 18 PetscCall(PetscOptionsFList("-mat_seqaij_type","Matrix SeqAIJ type","MatSeqAIJSetType",MatSeqAIJList,"seqaij",type,256,&flg)); 19 if (flg) PetscCall(MatSeqAIJSetType(A,type)); 20 PetscOptionsEnd(); 21 PetscFunctionReturn(0); 22 } 23 24 PetscErrorCode MatGetColumnReductions_SeqAIJ(Mat A,PetscInt type,PetscReal *reductions) 25 { 26 PetscInt i,m,n; 27 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 28 29 PetscFunctionBegin; 30 PetscCall(MatGetSize(A,&m,&n)); 31 PetscCall(PetscArrayzero(reductions,n)); 32 if (type == NORM_2) { 33 for (i=0; i<aij->i[m]; i++) { 34 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 35 } 36 } else if (type == NORM_1) { 37 for (i=0; i<aij->i[m]; i++) { 38 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]); 39 } 40 } else if (type == NORM_INFINITY) { 41 for (i=0; i<aij->i[m]; i++) { 42 reductions[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),reductions[aij->j[i]]); 43 } 44 } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) { 45 for (i=0; i<aij->i[m]; i++) { 46 reductions[aij->j[i]] += PetscRealPart(aij->a[i]); 47 } 48 } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) { 49 for (i=0; i<aij->i[m]; i++) { 50 reductions[aij->j[i]] += PetscImaginaryPart(aij->a[i]); 51 } 52 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown reduction type"); 53 54 if (type == NORM_2) { 55 for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]); 56 } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) { 57 for (i=0; i<n; i++) reductions[i] /= m; 58 } 59 PetscFunctionReturn(0); 60 } 61 62 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is) 63 { 64 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 65 PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs; 66 const PetscInt *jj = a->j,*ii = a->i; 67 PetscInt *rows; 68 69 PetscFunctionBegin; 70 for (i=0; i<m; i++) { 71 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 72 cnt++; 73 } 74 } 75 PetscCall(PetscMalloc1(cnt,&rows)); 76 cnt = 0; 77 for (i=0; i<m; i++) { 78 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 79 rows[cnt] = i; 80 cnt++; 81 } 82 } 83 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is)); 84 PetscFunctionReturn(0); 85 } 86 87 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 88 { 89 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 90 const MatScalar *aa; 91 PetscInt i,m=A->rmap->n,cnt = 0; 92 const PetscInt *ii = a->i,*jj = a->j,*diag; 93 PetscInt *rows; 94 95 PetscFunctionBegin; 96 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 97 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 98 diag = a->diag; 99 for (i=0; i<m; i++) { 100 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 101 cnt++; 102 } 103 } 104 PetscCall(PetscMalloc1(cnt,&rows)); 105 cnt = 0; 106 for (i=0; i<m; i++) { 107 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 108 rows[cnt++] = i; 109 } 110 } 111 *nrows = cnt; 112 *zrows = rows; 113 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 118 { 119 PetscInt nrows,*rows; 120 121 PetscFunctionBegin; 122 *zrows = NULL; 123 PetscCall(MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows)); 124 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows)); 125 PetscFunctionReturn(0); 126 } 127 128 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 129 { 130 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 131 const MatScalar *aa; 132 PetscInt m=A->rmap->n,cnt = 0; 133 const PetscInt *ii; 134 PetscInt n,i,j,*rows; 135 136 PetscFunctionBegin; 137 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 138 *keptrows = NULL; 139 ii = a->i; 140 for (i=0; i<m; i++) { 141 n = ii[i+1] - ii[i]; 142 if (!n) { 143 cnt++; 144 goto ok1; 145 } 146 for (j=ii[i]; j<ii[i+1]; j++) { 147 if (aa[j] != 0.0) goto ok1; 148 } 149 cnt++; 150 ok1:; 151 } 152 if (!cnt) { 153 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 154 PetscFunctionReturn(0); 155 } 156 PetscCall(PetscMalloc1(A->rmap->n-cnt,&rows)); 157 cnt = 0; 158 for (i=0; i<m; i++) { 159 n = ii[i+1] - ii[i]; 160 if (!n) continue; 161 for (j=ii[i]; j<ii[i+1]; j++) { 162 if (aa[j] != 0.0) { 163 rows[cnt++] = i; 164 break; 165 } 166 } 167 } 168 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 169 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows)); 170 PetscFunctionReturn(0); 171 } 172 173 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 174 { 175 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 176 PetscInt i,m = Y->rmap->n; 177 const PetscInt *diag; 178 MatScalar *aa; 179 const PetscScalar *v; 180 PetscBool missing; 181 182 PetscFunctionBegin; 183 if (Y->assembled) { 184 PetscCall(MatMissingDiagonal_SeqAIJ(Y,&missing,NULL)); 185 if (!missing) { 186 diag = aij->diag; 187 PetscCall(VecGetArrayRead(D,&v)); 188 PetscCall(MatSeqAIJGetArray(Y,&aa)); 189 if (is == INSERT_VALUES) { 190 for (i=0; i<m; i++) { 191 aa[diag[i]] = v[i]; 192 } 193 } else { 194 for (i=0; i<m; i++) { 195 aa[diag[i]] += v[i]; 196 } 197 } 198 PetscCall(MatSeqAIJRestoreArray(Y,&aa)); 199 PetscCall(VecRestoreArrayRead(D,&v)); 200 PetscFunctionReturn(0); 201 } 202 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 203 } 204 PetscCall(MatDiagonalSet_Default(Y,D,is)); 205 PetscFunctionReturn(0); 206 } 207 208 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 209 { 210 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 211 PetscInt i,ishift; 212 213 PetscFunctionBegin; 214 *m = A->rmap->n; 215 if (!ia) PetscFunctionReturn(0); 216 ishift = 0; 217 if (symmetric && !A->structurally_symmetric) { 218 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja)); 219 } else if (oshift == 1) { 220 PetscInt *tia; 221 PetscInt nz = a->i[A->rmap->n]; 222 /* malloc space and add 1 to i and j indices */ 223 PetscCall(PetscMalloc1(A->rmap->n+1,&tia)); 224 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 225 *ia = tia; 226 if (ja) { 227 PetscInt *tja; 228 PetscCall(PetscMalloc1(nz+1,&tja)); 229 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 230 *ja = tja; 231 } 232 } else { 233 *ia = a->i; 234 if (ja) *ja = a->j; 235 } 236 PetscFunctionReturn(0); 237 } 238 239 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 240 { 241 PetscFunctionBegin; 242 if (!ia) PetscFunctionReturn(0); 243 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 244 PetscCall(PetscFree(*ia)); 245 if (ja) PetscCall(PetscFree(*ja)); 246 } 247 PetscFunctionReturn(0); 248 } 249 250 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 251 { 252 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 253 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 254 PetscInt nz = a->i[m],row,*jj,mr,col; 255 256 PetscFunctionBegin; 257 *nn = n; 258 if (!ia) PetscFunctionReturn(0); 259 if (symmetric) { 260 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,0,oshift,(PetscInt**)ia,(PetscInt**)ja)); 261 } else { 262 PetscCall(PetscCalloc1(n,&collengths)); 263 PetscCall(PetscMalloc1(n+1,&cia)); 264 PetscCall(PetscMalloc1(nz,&cja)); 265 jj = a->j; 266 for (i=0; i<nz; i++) { 267 collengths[jj[i]]++; 268 } 269 cia[0] = oshift; 270 for (i=0; i<n; i++) { 271 cia[i+1] = cia[i] + collengths[i]; 272 } 273 PetscCall(PetscArrayzero(collengths,n)); 274 jj = a->j; 275 for (row=0; row<m; row++) { 276 mr = a->i[row+1] - a->i[row]; 277 for (i=0; i<mr; i++) { 278 col = *jj++; 279 280 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 281 } 282 } 283 PetscCall(PetscFree(collengths)); 284 *ia = cia; *ja = cja; 285 } 286 PetscFunctionReturn(0); 287 } 288 289 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 290 { 291 PetscFunctionBegin; 292 if (!ia) PetscFunctionReturn(0); 293 294 PetscCall(PetscFree(*ia)); 295 PetscCall(PetscFree(*ja)); 296 PetscFunctionReturn(0); 297 } 298 299 /* 300 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 301 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 302 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 303 */ 304 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 305 { 306 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 307 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 308 PetscInt nz = a->i[m],row,mr,col,tmp; 309 PetscInt *cspidx; 310 const PetscInt *jj; 311 312 PetscFunctionBegin; 313 *nn = n; 314 if (!ia) PetscFunctionReturn(0); 315 316 PetscCall(PetscCalloc1(n,&collengths)); 317 PetscCall(PetscMalloc1(n+1,&cia)); 318 PetscCall(PetscMalloc1(nz,&cja)); 319 PetscCall(PetscMalloc1(nz,&cspidx)); 320 jj = a->j; 321 for (i=0; i<nz; i++) { 322 collengths[jj[i]]++; 323 } 324 cia[0] = oshift; 325 for (i=0; i<n; i++) { 326 cia[i+1] = cia[i] + collengths[i]; 327 } 328 PetscCall(PetscArrayzero(collengths,n)); 329 jj = a->j; 330 for (row=0; row<m; row++) { 331 mr = a->i[row+1] - a->i[row]; 332 for (i=0; i<mr; i++) { 333 col = *jj++; 334 tmp = cia[col] + collengths[col]++ - oshift; 335 cspidx[tmp] = a->i[row] + i; /* index of a->j */ 336 cja[tmp] = row + oshift; 337 } 338 } 339 PetscCall(PetscFree(collengths)); 340 *ia = cia; 341 *ja = cja; 342 *spidx = cspidx; 343 PetscFunctionReturn(0); 344 } 345 346 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 347 { 348 PetscFunctionBegin; 349 PetscCall(MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done)); 350 PetscCall(PetscFree(*spidx)); 351 PetscFunctionReturn(0); 352 } 353 354 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 355 { 356 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 357 PetscInt *ai = a->i; 358 PetscScalar *aa; 359 360 PetscFunctionBegin; 361 PetscCall(MatSeqAIJGetArray(A,&aa)); 362 PetscCall(PetscArraycpy(aa+ai[row],v,ai[row+1]-ai[row])); 363 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 364 PetscFunctionReturn(0); 365 } 366 367 /* 368 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 369 370 - a single row of values is set with each call 371 - no row or column indices are negative or (in error) larger than the number of rows or columns 372 - the values are always added to the matrix, not set 373 - no new locations are introduced in the nonzero structure of the matrix 374 375 This does NOT assume the global column indices are sorted 376 377 */ 378 379 #include <petsc/private/isimpl.h> 380 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 381 { 382 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 383 PetscInt low,high,t,row,nrow,i,col,l; 384 const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j; 385 PetscInt lastcol = -1; 386 MatScalar *ap,value,*aa; 387 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 388 389 PetscFunctionBegin; 390 PetscCall(MatSeqAIJGetArray(A,&aa)); 391 row = ridx[im[0]]; 392 rp = aj + ai[row]; 393 ap = aa + ai[row]; 394 nrow = ailen[row]; 395 low = 0; 396 high = nrow; 397 for (l=0; l<n; l++) { /* loop over added columns */ 398 col = cidx[in[l]]; 399 value = v[l]; 400 401 if (col <= lastcol) low = 0; 402 else high = nrow; 403 lastcol = col; 404 while (high-low > 5) { 405 t = (low+high)/2; 406 if (rp[t] > col) high = t; 407 else low = t; 408 } 409 for (i=low; i<high; i++) { 410 if (rp[i] == col) { 411 ap[i] += value; 412 low = i + 1; 413 break; 414 } 415 } 416 } 417 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 418 return 0; 419 } 420 421 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 422 { 423 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 424 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 425 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 426 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 427 MatScalar *ap=NULL,value=0.0,*aa; 428 PetscBool ignorezeroentries = a->ignorezeroentries; 429 PetscBool roworiented = a->roworiented; 430 431 PetscFunctionBegin; 432 PetscCall(MatSeqAIJGetArray(A,&aa)); 433 for (k=0; k<m; k++) { /* loop over added rows */ 434 row = im[k]; 435 if (row < 0) continue; 436 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 437 rp = aj + ai[row]; 438 if (!A->structure_only) ap = aa + ai[row]; 439 rmax = imax[row]; nrow = ailen[row]; 440 low = 0; 441 high = nrow; 442 for (l=0; l<n; l++) { /* loop over added columns */ 443 if (in[l] < 0) continue; 444 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 445 col = in[l]; 446 if (v && !A->structure_only) value = roworiented ? v[l + k*n] : v[k + l*m]; 447 if (!A->structure_only && value == 0.0 && ignorezeroentries && is == ADD_VALUES && row != col) continue; 448 449 if (col <= lastcol) low = 0; 450 else high = nrow; 451 lastcol = col; 452 while (high-low > 5) { 453 t = (low+high)/2; 454 if (rp[t] > col) high = t; 455 else low = t; 456 } 457 for (i=low; i<high; i++) { 458 if (rp[i] > col) break; 459 if (rp[i] == col) { 460 if (!A->structure_only) { 461 if (is == ADD_VALUES) { 462 ap[i] += value; 463 (void)PetscLogFlops(1.0); 464 } 465 else ap[i] = value; 466 } 467 low = i + 1; 468 goto noinsert; 469 } 470 } 471 if (value == 0.0 && ignorezeroentries && row != col) goto noinsert; 472 if (nonew == 1) goto noinsert; 473 PetscCheck(nonew != -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%" PetscInt_FMT ",%" PetscInt_FMT ") in the matrix",row,col); 474 if (A->structure_only) { 475 MatSeqXAIJReallocateAIJ_structure_only(A,A->rmap->n,1,nrow,row,col,rmax,ai,aj,rp,imax,nonew,MatScalar); 476 } else { 477 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 478 } 479 N = nrow++ - 1; a->nz++; high++; 480 /* shift up all the later entries in this row */ 481 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 482 rp[i] = col; 483 if (!A->structure_only) { 484 PetscCall(PetscArraymove(ap+i+1,ap+i,N-i+1)); 485 ap[i] = value; 486 } 487 low = i + 1; 488 A->nonzerostate++; 489 noinsert:; 490 } 491 ailen[row] = nrow; 492 } 493 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 494 PetscFunctionReturn(0); 495 } 496 497 PetscErrorCode MatSetValues_SeqAIJ_SortedFullNoPreallocation(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 498 { 499 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 500 PetscInt *rp,k,row; 501 PetscInt *ai = a->i; 502 PetscInt *aj = a->j; 503 MatScalar *aa,*ap; 504 505 PetscFunctionBegin; 506 PetscCheck(!A->was_assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot call on assembled matrix."); 507 PetscCheck(m*n+a->nz <= a->maxnz,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of entries in matrix will be larger than maximum nonzeros allocated for %" PetscInt_FMT " in MatSeqAIJSetTotalPreallocation()",a->maxnz); 508 509 PetscCall(MatSeqAIJGetArray(A,&aa)); 510 for (k=0; k<m; k++) { /* loop over added rows */ 511 row = im[k]; 512 rp = aj + ai[row]; 513 ap = aa + ai[row]; 514 515 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 516 if (!A->structure_only) { 517 if (v) { 518 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 519 v += n; 520 } else { 521 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 522 } 523 } 524 a->ilen[row] = n; 525 a->imax[row] = n; 526 a->i[row+1] = a->i[row]+n; 527 a->nz += n; 528 } 529 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 530 PetscFunctionReturn(0); 531 } 532 533 /*@ 534 MatSeqAIJSetTotalPreallocation - Sets an upper bound on the total number of expected nonzeros in the matrix. 535 536 Input Parameters: 537 + A - the SeqAIJ matrix 538 - nztotal - bound on the number of nonzeros 539 540 Level: advanced 541 542 Notes: 543 This can be called if you will be provided the matrix row by row (from row zero) with sorted column indices for each row. 544 Simply call MatSetValues() after this call to provide the matrix entries in the usual manner. This matrix may be used 545 as always with multiple matrix assemblies. 546 547 .seealso: MatSetOption(), MAT_SORTED_FULL, MatSetValues(), MatSeqAIJSetPreallocation() 548 @*/ 549 550 PetscErrorCode MatSeqAIJSetTotalPreallocation(Mat A,PetscInt nztotal) 551 { 552 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 553 554 PetscFunctionBegin; 555 PetscCall(PetscLayoutSetUp(A->rmap)); 556 PetscCall(PetscLayoutSetUp(A->cmap)); 557 a->maxnz = nztotal; 558 if (!a->imax) { 559 PetscCall(PetscMalloc1(A->rmap->n,&a->imax)); 560 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 561 } 562 if (!a->ilen) { 563 PetscCall(PetscMalloc1(A->rmap->n,&a->ilen)); 564 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 565 } else { 566 PetscCall(PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt))); 567 } 568 569 /* allocate the matrix space */ 570 if (A->structure_only) { 571 PetscCall(PetscMalloc1(nztotal,&a->j)); 572 PetscCall(PetscMalloc1(A->rmap->n+1,&a->i)); 573 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*sizeof(PetscInt))); 574 } else { 575 PetscCall(PetscMalloc3(nztotal,&a->a,nztotal,&a->j,A->rmap->n+1,&a->i)); 576 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*(sizeof(PetscScalar)+sizeof(PetscInt)))); 577 } 578 a->i[0] = 0; 579 if (A->structure_only) { 580 a->singlemalloc = PETSC_FALSE; 581 a->free_a = PETSC_FALSE; 582 } else { 583 a->singlemalloc = PETSC_TRUE; 584 a->free_a = PETSC_TRUE; 585 } 586 a->free_ij = PETSC_TRUE; 587 A->ops->setvalues = MatSetValues_SeqAIJ_SortedFullNoPreallocation; 588 A->preallocated = PETSC_TRUE; 589 PetscFunctionReturn(0); 590 } 591 592 PetscErrorCode MatSetValues_SeqAIJ_SortedFull(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 593 { 594 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 595 PetscInt *rp,k,row; 596 PetscInt *ai = a->i,*ailen = a->ilen; 597 PetscInt *aj = a->j; 598 MatScalar *aa,*ap; 599 600 PetscFunctionBegin; 601 PetscCall(MatSeqAIJGetArray(A,&aa)); 602 for (k=0; k<m; k++) { /* loop over added rows */ 603 row = im[k]; 604 PetscCheck(n <= a->imax[row],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Preallocation for row %" PetscInt_FMT " does not match number of columns provided",n); 605 rp = aj + ai[row]; 606 ap = aa + ai[row]; 607 if (!A->was_assembled) { 608 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 609 } 610 if (!A->structure_only) { 611 if (v) { 612 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 613 v += n; 614 } else { 615 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 616 } 617 } 618 ailen[row] = n; 619 a->nz += n; 620 } 621 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 622 PetscFunctionReturn(0); 623 } 624 625 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 626 { 627 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 628 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 629 PetscInt *ai = a->i,*ailen = a->ilen; 630 MatScalar *ap,*aa; 631 632 PetscFunctionBegin; 633 PetscCall(MatSeqAIJGetArray(A,&aa)); 634 for (k=0; k<m; k++) { /* loop over rows */ 635 row = im[k]; 636 if (row < 0) {v += n; continue;} /* negative row */ 637 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 638 rp = aj + ai[row]; ap = aa + ai[row]; 639 nrow = ailen[row]; 640 for (l=0; l<n; l++) { /* loop over columns */ 641 if (in[l] < 0) {v++; continue;} /* negative column */ 642 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 643 col = in[l]; 644 high = nrow; low = 0; /* assume unsorted */ 645 while (high-low > 5) { 646 t = (low+high)/2; 647 if (rp[t] > col) high = t; 648 else low = t; 649 } 650 for (i=low; i<high; i++) { 651 if (rp[i] > col) break; 652 if (rp[i] == col) { 653 *v++ = ap[i]; 654 goto finished; 655 } 656 } 657 *v++ = 0.0; 658 finished:; 659 } 660 } 661 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 662 PetscFunctionReturn(0); 663 } 664 665 PetscErrorCode MatView_SeqAIJ_Binary(Mat mat,PetscViewer viewer) 666 { 667 Mat_SeqAIJ *A = (Mat_SeqAIJ*)mat->data; 668 const PetscScalar *av; 669 PetscInt header[4],M,N,m,nz,i; 670 PetscInt *rowlens; 671 672 PetscFunctionBegin; 673 PetscCall(PetscViewerSetUp(viewer)); 674 675 M = mat->rmap->N; 676 N = mat->cmap->N; 677 m = mat->rmap->n; 678 nz = A->nz; 679 680 /* write matrix header */ 681 header[0] = MAT_FILE_CLASSID; 682 header[1] = M; header[2] = N; header[3] = nz; 683 PetscCall(PetscViewerBinaryWrite(viewer,header,4,PETSC_INT)); 684 685 /* fill in and store row lengths */ 686 PetscCall(PetscMalloc1(m,&rowlens)); 687 for (i=0; i<m; i++) rowlens[i] = A->i[i+1] - A->i[i]; 688 PetscCall(PetscViewerBinaryWrite(viewer,rowlens,m,PETSC_INT)); 689 PetscCall(PetscFree(rowlens)); 690 /* store column indices */ 691 PetscCall(PetscViewerBinaryWrite(viewer,A->j,nz,PETSC_INT)); 692 /* store nonzero values */ 693 PetscCall(MatSeqAIJGetArrayRead(mat,&av)); 694 PetscCall(PetscViewerBinaryWrite(viewer,av,nz,PETSC_SCALAR)); 695 PetscCall(MatSeqAIJRestoreArrayRead(mat,&av)); 696 697 /* write block size option to the viewer's .info file */ 698 PetscCall(MatView_Binary_BlockSizes(mat,viewer)); 699 PetscFunctionReturn(0); 700 } 701 702 static PetscErrorCode MatView_SeqAIJ_ASCII_structonly(Mat A,PetscViewer viewer) 703 { 704 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 705 PetscInt i,k,m=A->rmap->N; 706 707 PetscFunctionBegin; 708 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 709 for (i=0; i<m; i++) { 710 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 711 for (k=a->i[i]; k<a->i[i+1]; k++) { 712 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ") ",a->j[k])); 713 } 714 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 715 } 716 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 717 PetscFunctionReturn(0); 718 } 719 720 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 721 722 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 723 { 724 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 725 const PetscScalar *av; 726 PetscInt i,j,m = A->rmap->n; 727 const char *name; 728 PetscViewerFormat format; 729 730 PetscFunctionBegin; 731 if (A->structure_only) { 732 PetscCall(MatView_SeqAIJ_ASCII_structonly(A,viewer)); 733 PetscFunctionReturn(0); 734 } 735 736 PetscCall(PetscViewerGetFormat(viewer,&format)); 737 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 738 739 /* trigger copy to CPU if needed */ 740 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 741 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 742 if (format == PETSC_VIEWER_ASCII_MATLAB) { 743 PetscInt nofinalvalue = 0; 744 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) { 745 /* Need a dummy value to ensure the dimension of the matrix. */ 746 nofinalvalue = 1; 747 } 748 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 749 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n",m,A->cmap->n)); 750 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %" PetscInt_FMT " \n",a->nz)); 751 #if defined(PETSC_USE_COMPLEX) 752 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",4);\n",a->nz+nofinalvalue)); 753 #else 754 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",a->nz+nofinalvalue)); 755 #endif 756 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n")); 757 758 for (i=0; i<m; i++) { 759 for (j=a->i[i]; j<a->i[i+1]; j++) { 760 #if defined(PETSC_USE_COMPLEX) 761 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",i+1,a->j[j]+1,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 762 #else 763 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,a->j[j]+1,(double)a->a[j])); 764 #endif 765 } 766 } 767 if (nofinalvalue) { 768 #if defined(PETSC_USE_COMPLEX) 769 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",m,A->cmap->n,0.,0.)); 770 #else 771 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",m,A->cmap->n,0.0)); 772 #endif 773 } 774 PetscCall(PetscObjectGetName((PetscObject)A,&name)); 775 PetscCall(PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name)); 776 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 777 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 778 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 779 for (i=0; i<m; i++) { 780 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 781 for (j=a->i[i]; j<a->i[i+1]; j++) { 782 #if defined(PETSC_USE_COMPLEX) 783 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 784 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 785 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 786 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 787 } else if (PetscRealPart(a->a[j]) != 0.0) { 788 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 789 } 790 #else 791 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 792 #endif 793 } 794 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 795 } 796 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 797 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 798 PetscInt nzd=0,fshift=1,*sptr; 799 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 800 PetscCall(PetscMalloc1(m+1,&sptr)); 801 for (i=0; i<m; i++) { 802 sptr[i] = nzd+1; 803 for (j=a->i[i]; j<a->i[i+1]; j++) { 804 if (a->j[j] >= i) { 805 #if defined(PETSC_USE_COMPLEX) 806 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 807 #else 808 if (a->a[j] != 0.0) nzd++; 809 #endif 810 } 811 } 812 } 813 sptr[m] = nzd+1; 814 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n\n",m,nzd)); 815 for (i=0; i<m+1; i+=6) { 816 if (i+4<m) { 817 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5])); 818 } else if (i+3<m) { 819 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4])); 820 } else if (i+2<m) { 821 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3])); 822 } else if (i+1<m) { 823 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2])); 824 } else if (i<m) { 825 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1])); 826 } else { 827 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT "\n",sptr[i])); 828 } 829 } 830 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 831 PetscCall(PetscFree(sptr)); 832 for (i=0; i<m; i++) { 833 for (j=a->i[i]; j<a->i[i+1]; j++) { 834 if (a->j[j] >= i) PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " ",a->j[j]+fshift)); 835 } 836 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 837 } 838 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 839 for (i=0; i<m; i++) { 840 for (j=a->i[i]; j<a->i[i+1]; j++) { 841 if (a->j[j] >= i) { 842 #if defined(PETSC_USE_COMPLEX) 843 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 844 PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 845 } 846 #else 847 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j])); 848 #endif 849 } 850 } 851 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 852 } 853 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 854 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 855 PetscInt cnt = 0,jcnt; 856 PetscScalar value; 857 #if defined(PETSC_USE_COMPLEX) 858 PetscBool realonly = PETSC_TRUE; 859 860 for (i=0; i<a->i[m]; i++) { 861 if (PetscImaginaryPart(a->a[i]) != 0.0) { 862 realonly = PETSC_FALSE; 863 break; 864 } 865 } 866 #endif 867 868 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 869 for (i=0; i<m; i++) { 870 jcnt = 0; 871 for (j=0; j<A->cmap->n; j++) { 872 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 873 value = a->a[cnt++]; 874 jcnt++; 875 } else { 876 value = 0.0; 877 } 878 #if defined(PETSC_USE_COMPLEX) 879 if (realonly) { 880 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value))); 881 } else { 882 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value))); 883 } 884 #else 885 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value)); 886 #endif 887 } 888 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 889 } 890 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 891 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 892 PetscInt fshift=1; 893 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 894 #if defined(PETSC_USE_COMPLEX) 895 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n")); 896 #else 897 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n")); 898 #endif 899 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", m, A->cmap->n, a->nz)); 900 for (i=0; i<m; i++) { 901 for (j=a->i[i]; j<a->i[i+1]; j++) { 902 #if defined(PETSC_USE_COMPLEX) 903 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g %g\n", i+fshift,a->j[j]+fshift,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 904 #else 905 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j])); 906 #endif 907 } 908 } 909 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 910 } else { 911 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 912 if (A->factortype) { 913 for (i=0; i<m; i++) { 914 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 915 /* L part */ 916 for (j=a->i[i]; j<a->i[i+1]; j++) { 917 #if defined(PETSC_USE_COMPLEX) 918 if (PetscImaginaryPart(a->a[j]) > 0.0) { 919 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 920 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 921 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 922 } else { 923 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 924 } 925 #else 926 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 927 #endif 928 } 929 /* diagonal */ 930 j = a->diag[i]; 931 #if defined(PETSC_USE_COMPLEX) 932 if (PetscImaginaryPart(a->a[j]) > 0.0) { 933 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)PetscImaginaryPart(1.0/a->a[j]))); 934 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 935 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)(-PetscImaginaryPart(1.0/a->a[j])))); 936 } else { 937 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]))); 938 } 939 #else 940 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)(1.0/a->a[j]))); 941 #endif 942 943 /* U part */ 944 for (j=a->diag[i+1]+1; j<a->diag[i]; j++) { 945 #if defined(PETSC_USE_COMPLEX) 946 if (PetscImaginaryPart(a->a[j]) > 0.0) { 947 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 948 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 949 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 950 } else { 951 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 952 } 953 #else 954 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 955 #endif 956 } 957 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 958 } 959 } else { 960 for (i=0; i<m; i++) { 961 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 962 for (j=a->i[i]; j<a->i[i+1]; j++) { 963 #if defined(PETSC_USE_COMPLEX) 964 if (PetscImaginaryPart(a->a[j]) > 0.0) { 965 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 966 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 967 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 968 } else { 969 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 970 } 971 #else 972 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 973 #endif 974 } 975 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 976 } 977 } 978 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 979 } 980 PetscCall(PetscViewerFlush(viewer)); 981 PetscFunctionReturn(0); 982 } 983 984 #include <petscdraw.h> 985 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 986 { 987 Mat A = (Mat) Aa; 988 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 989 PetscInt i,j,m = A->rmap->n; 990 int color; 991 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 992 PetscViewer viewer; 993 PetscViewerFormat format; 994 const PetscScalar *aa; 995 996 PetscFunctionBegin; 997 PetscCall(PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer)); 998 PetscCall(PetscViewerGetFormat(viewer,&format)); 999 PetscCall(PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr)); 1000 1001 /* loop over matrix elements drawing boxes */ 1002 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1003 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1004 PetscDrawCollectiveBegin(draw); 1005 /* Blue for negative, Cyan for zero and Red for positive */ 1006 color = PETSC_DRAW_BLUE; 1007 for (i=0; i<m; i++) { 1008 y_l = m - i - 1.0; y_r = y_l + 1.0; 1009 for (j=a->i[i]; j<a->i[i+1]; j++) { 1010 x_l = a->j[j]; x_r = x_l + 1.0; 1011 if (PetscRealPart(aa[j]) >= 0.) continue; 1012 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1013 } 1014 } 1015 color = PETSC_DRAW_CYAN; 1016 for (i=0; i<m; i++) { 1017 y_l = m - i - 1.0; y_r = y_l + 1.0; 1018 for (j=a->i[i]; j<a->i[i+1]; j++) { 1019 x_l = a->j[j]; x_r = x_l + 1.0; 1020 if (aa[j] != 0.) continue; 1021 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1022 } 1023 } 1024 color = PETSC_DRAW_RED; 1025 for (i=0; i<m; i++) { 1026 y_l = m - i - 1.0; y_r = y_l + 1.0; 1027 for (j=a->i[i]; j<a->i[i+1]; j++) { 1028 x_l = a->j[j]; x_r = x_l + 1.0; 1029 if (PetscRealPart(aa[j]) <= 0.) continue; 1030 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1031 } 1032 } 1033 PetscDrawCollectiveEnd(draw); 1034 } else { 1035 /* use contour shading to indicate magnitude of values */ 1036 /* first determine max of all nonzero values */ 1037 PetscReal minv = 0.0, maxv = 0.0; 1038 PetscInt nz = a->nz, count = 0; 1039 PetscDraw popup; 1040 1041 for (i=0; i<nz; i++) { 1042 if (PetscAbsScalar(aa[i]) > maxv) maxv = PetscAbsScalar(aa[i]); 1043 } 1044 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1045 PetscCall(PetscDrawGetPopup(draw,&popup)); 1046 PetscCall(PetscDrawScalePopup(popup,minv,maxv)); 1047 1048 PetscDrawCollectiveBegin(draw); 1049 for (i=0; i<m; i++) { 1050 y_l = m - i - 1.0; 1051 y_r = y_l + 1.0; 1052 for (j=a->i[i]; j<a->i[i+1]; j++) { 1053 x_l = a->j[j]; 1054 x_r = x_l + 1.0; 1055 color = PetscDrawRealToColor(PetscAbsScalar(aa[count]),minv,maxv); 1056 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1057 count++; 1058 } 1059 } 1060 PetscDrawCollectiveEnd(draw); 1061 } 1062 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1063 PetscFunctionReturn(0); 1064 } 1065 1066 #include <petscdraw.h> 1067 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 1068 { 1069 PetscDraw draw; 1070 PetscReal xr,yr,xl,yl,h,w; 1071 PetscBool isnull; 1072 1073 PetscFunctionBegin; 1074 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 1075 PetscCall(PetscDrawIsNull(draw,&isnull)); 1076 if (isnull) PetscFunctionReturn(0); 1077 1078 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1079 xr += w; yr += h; xl = -w; yl = -h; 1080 PetscCall(PetscDrawSetCoordinates(draw,xl,yl,xr,yr)); 1081 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer)); 1082 PetscCall(PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A)); 1083 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL)); 1084 PetscCall(PetscDrawSave(draw)); 1085 PetscFunctionReturn(0); 1086 } 1087 1088 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 1089 { 1090 PetscBool iascii,isbinary,isdraw; 1091 1092 PetscFunctionBegin; 1093 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1094 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 1095 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 1096 if (iascii) { 1097 PetscCall(MatView_SeqAIJ_ASCII(A,viewer)); 1098 } else if (isbinary) { 1099 PetscCall(MatView_SeqAIJ_Binary(A,viewer)); 1100 } else if (isdraw) { 1101 PetscCall(MatView_SeqAIJ_Draw(A,viewer)); 1102 } 1103 PetscCall(MatView_SeqAIJ_Inode(A,viewer)); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 1108 { 1109 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1110 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 1111 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 1112 MatScalar *aa = a->a,*ap; 1113 PetscReal ratio = 0.6; 1114 1115 PetscFunctionBegin; 1116 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 1117 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1118 if (A->was_assembled && A->ass_nonzerostate == A->nonzerostate) { 1119 /* we need to respect users asking to use or not the inodes routine in between matrix assemblies */ 1120 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1121 PetscFunctionReturn(0); 1122 } 1123 1124 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 1125 for (i=1; i<m; i++) { 1126 /* move each row back by the amount of empty slots (fshift) before it*/ 1127 fshift += imax[i-1] - ailen[i-1]; 1128 rmax = PetscMax(rmax,ailen[i]); 1129 if (fshift) { 1130 ip = aj + ai[i]; 1131 ap = aa + ai[i]; 1132 N = ailen[i]; 1133 PetscCall(PetscArraymove(ip-fshift,ip,N)); 1134 if (!A->structure_only) { 1135 PetscCall(PetscArraymove(ap-fshift,ap,N)); 1136 } 1137 } 1138 ai[i] = ai[i-1] + ailen[i-1]; 1139 } 1140 if (m) { 1141 fshift += imax[m-1] - ailen[m-1]; 1142 ai[m] = ai[m-1] + ailen[m-1]; 1143 } 1144 1145 /* reset ilen and imax for each row */ 1146 a->nonzerorowcnt = 0; 1147 if (A->structure_only) { 1148 PetscCall(PetscFree(a->imax)); 1149 PetscCall(PetscFree(a->ilen)); 1150 } else { /* !A->structure_only */ 1151 for (i=0; i<m; i++) { 1152 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1153 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 1154 } 1155 } 1156 a->nz = ai[m]; 1157 PetscCheckFalse(fshift && a->nounused == -1,PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, fshift); 1158 1159 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1160 PetscCall(PetscInfo(A,"Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded,%" PetscInt_FMT " used\n",m,A->cmap->n,fshift,a->nz)); 1161 PetscCall(PetscInfo(A,"Number of mallocs during MatSetValues() is %" PetscInt_FMT "\n",a->reallocs)); 1162 PetscCall(PetscInfo(A,"Maximum nonzeros in any row is %" PetscInt_FMT "\n",rmax)); 1163 1164 A->info.mallocs += a->reallocs; 1165 a->reallocs = 0; 1166 A->info.nz_unneeded = (PetscReal)fshift; 1167 a->rmax = rmax; 1168 1169 if (!A->structure_only) { 1170 PetscCall(MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio)); 1171 } 1172 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1173 PetscFunctionReturn(0); 1174 } 1175 1176 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1177 { 1178 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1179 PetscInt i,nz = a->nz; 1180 MatScalar *aa; 1181 1182 PetscFunctionBegin; 1183 PetscCall(MatSeqAIJGetArray(A,&aa)); 1184 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1185 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1186 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1187 PetscFunctionReturn(0); 1188 } 1189 1190 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1191 { 1192 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1193 PetscInt i,nz = a->nz; 1194 MatScalar *aa; 1195 1196 PetscFunctionBegin; 1197 PetscCall(MatSeqAIJGetArray(A,&aa)); 1198 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1199 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1200 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1201 PetscFunctionReturn(0); 1202 } 1203 1204 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1205 { 1206 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1207 MatScalar *aa; 1208 1209 PetscFunctionBegin; 1210 PetscCall(MatSeqAIJGetArrayWrite(A,&aa)); 1211 PetscCall(PetscArrayzero(aa,a->i[A->rmap->n])); 1212 PetscCall(MatSeqAIJRestoreArrayWrite(A,&aa)); 1213 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1214 PetscFunctionReturn(0); 1215 } 1216 1217 PETSC_INTERN PetscErrorCode MatResetPreallocationCOO_SeqAIJ(Mat A) 1218 { 1219 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1220 1221 PetscFunctionBegin; 1222 PetscCall(PetscFree(a->perm)); 1223 PetscCall(PetscFree(a->jmap)); 1224 PetscFunctionReturn(0); 1225 } 1226 1227 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1228 { 1229 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1230 1231 PetscFunctionBegin; 1232 #if defined(PETSC_USE_LOG) 1233 PetscLogObjectState((PetscObject)A,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT,A->rmap->n,A->cmap->n,a->nz); 1234 #endif 1235 PetscCall(MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i)); 1236 PetscCall(MatResetPreallocationCOO_SeqAIJ(A)); 1237 PetscCall(ISDestroy(&a->row)); 1238 PetscCall(ISDestroy(&a->col)); 1239 PetscCall(PetscFree(a->diag)); 1240 PetscCall(PetscFree(a->ibdiag)); 1241 PetscCall(PetscFree(a->imax)); 1242 PetscCall(PetscFree(a->ilen)); 1243 PetscCall(PetscFree(a->ipre)); 1244 PetscCall(PetscFree3(a->idiag,a->mdiag,a->ssor_work)); 1245 PetscCall(PetscFree(a->solve_work)); 1246 PetscCall(ISDestroy(&a->icol)); 1247 PetscCall(PetscFree(a->saved_values)); 1248 PetscCall(PetscFree2(a->compressedrow.i,a->compressedrow.rindex)); 1249 PetscCall(MatDestroy_SeqAIJ_Inode(A)); 1250 PetscCall(PetscFree(A->data)); 1251 1252 /* MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted may allocate this. 1253 That function is so heavily used (sometimes in an hidden way through multnumeric function pointers) 1254 that is hard to properly add this data to the MatProduct data. We free it here to avoid 1255 users reusing the matrix object with different data to incur in obscure segmentation faults 1256 due to different matrix sizes */ 1257 PetscCall(PetscObjectCompose((PetscObject)A,"__PETSc__ab_dense",NULL)); 1258 1259 PetscCall(PetscObjectChangeTypeName((PetscObject)A,NULL)); 1260 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL)); 1261 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL)); 1262 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL)); 1263 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL)); 1264 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL)); 1265 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL)); 1266 #if defined(PETSC_HAVE_CUDA) 1267 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcusparse_C",NULL)); 1268 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",NULL)); 1269 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",NULL)); 1270 #endif 1271 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 1272 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijkokkos_C",NULL)); 1273 #endif 1274 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcrl_C",NULL)); 1275 #if defined(PETSC_HAVE_ELEMENTAL) 1276 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL)); 1277 #endif 1278 #if defined(PETSC_HAVE_SCALAPACK) 1279 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_scalapack_C",NULL)); 1280 #endif 1281 #if defined(PETSC_HAVE_HYPRE) 1282 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_hypre_C",NULL)); 1283 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",NULL)); 1284 #endif 1285 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL)); 1286 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsell_C",NULL)); 1287 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_is_C",NULL)); 1288 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL)); 1289 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL)); 1290 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatResetPreallocation_C",NULL)); 1291 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL)); 1292 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL)); 1293 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_is_seqaij_C",NULL)); 1294 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqdense_seqaij_C",NULL)); 1295 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaij_C",NULL)); 1296 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJKron_C",NULL)); 1297 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL)); 1298 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL)); 1299 PetscFunctionReturn(0); 1300 } 1301 1302 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1303 { 1304 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1305 1306 PetscFunctionBegin; 1307 switch (op) { 1308 case MAT_ROW_ORIENTED: 1309 a->roworiented = flg; 1310 break; 1311 case MAT_KEEP_NONZERO_PATTERN: 1312 a->keepnonzeropattern = flg; 1313 break; 1314 case MAT_NEW_NONZERO_LOCATIONS: 1315 a->nonew = (flg ? 0 : 1); 1316 break; 1317 case MAT_NEW_NONZERO_LOCATION_ERR: 1318 a->nonew = (flg ? -1 : 0); 1319 break; 1320 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1321 a->nonew = (flg ? -2 : 0); 1322 break; 1323 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1324 a->nounused = (flg ? -1 : 0); 1325 break; 1326 case MAT_IGNORE_ZERO_ENTRIES: 1327 a->ignorezeroentries = flg; 1328 break; 1329 case MAT_SPD: 1330 case MAT_SYMMETRIC: 1331 case MAT_STRUCTURALLY_SYMMETRIC: 1332 case MAT_HERMITIAN: 1333 case MAT_SYMMETRY_ETERNAL: 1334 case MAT_STRUCTURE_ONLY: 1335 /* These options are handled directly by MatSetOption() */ 1336 break; 1337 case MAT_FORCE_DIAGONAL_ENTRIES: 1338 case MAT_IGNORE_OFF_PROC_ENTRIES: 1339 case MAT_USE_HASH_TABLE: 1340 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 1341 break; 1342 case MAT_USE_INODES: 1343 PetscCall(MatSetOption_SeqAIJ_Inode(A,MAT_USE_INODES,flg)); 1344 break; 1345 case MAT_SUBMAT_SINGLEIS: 1346 A->submat_singleis = flg; 1347 break; 1348 case MAT_SORTED_FULL: 1349 if (flg) A->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 1350 else A->ops->setvalues = MatSetValues_SeqAIJ; 1351 break; 1352 case MAT_FORM_EXPLICIT_TRANSPOSE: 1353 A->form_explicit_transpose = flg; 1354 break; 1355 default: 1356 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1357 } 1358 PetscFunctionReturn(0); 1359 } 1360 1361 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1362 { 1363 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1364 PetscInt i,j,n,*ai=a->i,*aj=a->j; 1365 PetscScalar *x; 1366 const PetscScalar *aa; 1367 1368 PetscFunctionBegin; 1369 PetscCall(VecGetLocalSize(v,&n)); 1370 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1371 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1372 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1373 PetscInt *diag=a->diag; 1374 PetscCall(VecGetArrayWrite(v,&x)); 1375 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1376 PetscCall(VecRestoreArrayWrite(v,&x)); 1377 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1378 PetscFunctionReturn(0); 1379 } 1380 1381 PetscCall(VecGetArrayWrite(v,&x)); 1382 for (i=0; i<n; i++) { 1383 x[i] = 0.0; 1384 for (j=ai[i]; j<ai[i+1]; j++) { 1385 if (aj[j] == i) { 1386 x[i] = aa[j]; 1387 break; 1388 } 1389 } 1390 } 1391 PetscCall(VecRestoreArrayWrite(v,&x)); 1392 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1393 PetscFunctionReturn(0); 1394 } 1395 1396 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1397 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1398 { 1399 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1400 const MatScalar *aa; 1401 PetscScalar *y; 1402 const PetscScalar *x; 1403 PetscInt m = A->rmap->n; 1404 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1405 const MatScalar *v; 1406 PetscScalar alpha; 1407 PetscInt n,i,j; 1408 const PetscInt *idx,*ii,*ridx=NULL; 1409 Mat_CompressedRow cprow = a->compressedrow; 1410 PetscBool usecprow = cprow.use; 1411 #endif 1412 1413 PetscFunctionBegin; 1414 if (zz != yy) PetscCall(VecCopy(zz,yy)); 1415 PetscCall(VecGetArrayRead(xx,&x)); 1416 PetscCall(VecGetArray(yy,&y)); 1417 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1418 1419 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1420 fortranmulttransposeaddaij_(&m,x,a->i,a->j,aa,y); 1421 #else 1422 if (usecprow) { 1423 m = cprow.nrows; 1424 ii = cprow.i; 1425 ridx = cprow.rindex; 1426 } else { 1427 ii = a->i; 1428 } 1429 for (i=0; i<m; i++) { 1430 idx = a->j + ii[i]; 1431 v = aa + ii[i]; 1432 n = ii[i+1] - ii[i]; 1433 if (usecprow) { 1434 alpha = x[ridx[i]]; 1435 } else { 1436 alpha = x[i]; 1437 } 1438 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1439 } 1440 #endif 1441 PetscCall(PetscLogFlops(2.0*a->nz)); 1442 PetscCall(VecRestoreArrayRead(xx,&x)); 1443 PetscCall(VecRestoreArray(yy,&y)); 1444 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1445 PetscFunctionReturn(0); 1446 } 1447 1448 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1449 { 1450 PetscFunctionBegin; 1451 PetscCall(VecSet(yy,0.0)); 1452 PetscCall(MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy)); 1453 PetscFunctionReturn(0); 1454 } 1455 1456 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1457 1458 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1459 { 1460 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1461 PetscScalar *y; 1462 const PetscScalar *x; 1463 const MatScalar *aa,*a_a; 1464 PetscInt m=A->rmap->n; 1465 const PetscInt *aj,*ii,*ridx=NULL; 1466 PetscInt n,i; 1467 PetscScalar sum; 1468 PetscBool usecprow=a->compressedrow.use; 1469 1470 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1471 #pragma disjoint(*x,*y,*aa) 1472 #endif 1473 1474 PetscFunctionBegin; 1475 if (a->inode.use && a->inode.checked) { 1476 PetscCall(MatMult_SeqAIJ_Inode(A,xx,yy)); 1477 PetscFunctionReturn(0); 1478 } 1479 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1480 PetscCall(VecGetArrayRead(xx,&x)); 1481 PetscCall(VecGetArray(yy,&y)); 1482 ii = a->i; 1483 if (usecprow) { /* use compressed row format */ 1484 PetscCall(PetscArrayzero(y,m)); 1485 m = a->compressedrow.nrows; 1486 ii = a->compressedrow.i; 1487 ridx = a->compressedrow.rindex; 1488 for (i=0; i<m; i++) { 1489 n = ii[i+1] - ii[i]; 1490 aj = a->j + ii[i]; 1491 aa = a_a + ii[i]; 1492 sum = 0.0; 1493 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1494 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1495 y[*ridx++] = sum; 1496 } 1497 } else { /* do not use compressed row format */ 1498 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1499 aj = a->j; 1500 aa = a_a; 1501 fortranmultaij_(&m,x,ii,aj,aa,y); 1502 #else 1503 for (i=0; i<m; i++) { 1504 n = ii[i+1] - ii[i]; 1505 aj = a->j + ii[i]; 1506 aa = a_a + ii[i]; 1507 sum = 0.0; 1508 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1509 y[i] = sum; 1510 } 1511 #endif 1512 } 1513 PetscCall(PetscLogFlops(2.0*a->nz - a->nonzerorowcnt)); 1514 PetscCall(VecRestoreArrayRead(xx,&x)); 1515 PetscCall(VecRestoreArray(yy,&y)); 1516 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1517 PetscFunctionReturn(0); 1518 } 1519 1520 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1521 { 1522 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1523 PetscScalar *y; 1524 const PetscScalar *x; 1525 const MatScalar *aa,*a_a; 1526 PetscInt m=A->rmap->n; 1527 const PetscInt *aj,*ii,*ridx=NULL; 1528 PetscInt n,i,nonzerorow=0; 1529 PetscScalar sum; 1530 PetscBool usecprow=a->compressedrow.use; 1531 1532 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1533 #pragma disjoint(*x,*y,*aa) 1534 #endif 1535 1536 PetscFunctionBegin; 1537 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1538 PetscCall(VecGetArrayRead(xx,&x)); 1539 PetscCall(VecGetArray(yy,&y)); 1540 if (usecprow) { /* use compressed row format */ 1541 m = a->compressedrow.nrows; 1542 ii = a->compressedrow.i; 1543 ridx = a->compressedrow.rindex; 1544 for (i=0; i<m; i++) { 1545 n = ii[i+1] - ii[i]; 1546 aj = a->j + ii[i]; 1547 aa = a_a + ii[i]; 1548 sum = 0.0; 1549 nonzerorow += (n>0); 1550 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1551 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1552 y[*ridx++] = sum; 1553 } 1554 } else { /* do not use compressed row format */ 1555 ii = a->i; 1556 for (i=0; i<m; i++) { 1557 n = ii[i+1] - ii[i]; 1558 aj = a->j + ii[i]; 1559 aa = a_a + ii[i]; 1560 sum = 0.0; 1561 nonzerorow += (n>0); 1562 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1563 y[i] = sum; 1564 } 1565 } 1566 PetscCall(PetscLogFlops(2.0*a->nz - nonzerorow)); 1567 PetscCall(VecRestoreArrayRead(xx,&x)); 1568 PetscCall(VecRestoreArray(yy,&y)); 1569 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1570 PetscFunctionReturn(0); 1571 } 1572 1573 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1574 { 1575 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1576 PetscScalar *y,*z; 1577 const PetscScalar *x; 1578 const MatScalar *aa,*a_a; 1579 PetscInt m = A->rmap->n,*aj,*ii; 1580 PetscInt n,i,*ridx=NULL; 1581 PetscScalar sum; 1582 PetscBool usecprow=a->compressedrow.use; 1583 1584 PetscFunctionBegin; 1585 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1586 PetscCall(VecGetArrayRead(xx,&x)); 1587 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1588 if (usecprow) { /* use compressed row format */ 1589 if (zz != yy) { 1590 PetscCall(PetscArraycpy(z,y,m)); 1591 } 1592 m = a->compressedrow.nrows; 1593 ii = a->compressedrow.i; 1594 ridx = a->compressedrow.rindex; 1595 for (i=0; i<m; i++) { 1596 n = ii[i+1] - ii[i]; 1597 aj = a->j + ii[i]; 1598 aa = a_a + ii[i]; 1599 sum = y[*ridx]; 1600 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1601 z[*ridx++] = sum; 1602 } 1603 } else { /* do not use compressed row format */ 1604 ii = a->i; 1605 for (i=0; i<m; i++) { 1606 n = ii[i+1] - ii[i]; 1607 aj = a->j + ii[i]; 1608 aa = a_a + ii[i]; 1609 sum = y[i]; 1610 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1611 z[i] = sum; 1612 } 1613 } 1614 PetscCall(PetscLogFlops(2.0*a->nz)); 1615 PetscCall(VecRestoreArrayRead(xx,&x)); 1616 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1617 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1618 PetscFunctionReturn(0); 1619 } 1620 1621 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1622 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1623 { 1624 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1625 PetscScalar *y,*z; 1626 const PetscScalar *x; 1627 const MatScalar *aa,*a_a; 1628 const PetscInt *aj,*ii,*ridx=NULL; 1629 PetscInt m = A->rmap->n,n,i; 1630 PetscScalar sum; 1631 PetscBool usecprow=a->compressedrow.use; 1632 1633 PetscFunctionBegin; 1634 if (a->inode.use && a->inode.checked) { 1635 PetscCall(MatMultAdd_SeqAIJ_Inode(A,xx,yy,zz)); 1636 PetscFunctionReturn(0); 1637 } 1638 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1639 PetscCall(VecGetArrayRead(xx,&x)); 1640 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1641 if (usecprow) { /* use compressed row format */ 1642 if (zz != yy) { 1643 PetscCall(PetscArraycpy(z,y,m)); 1644 } 1645 m = a->compressedrow.nrows; 1646 ii = a->compressedrow.i; 1647 ridx = a->compressedrow.rindex; 1648 for (i=0; i<m; i++) { 1649 n = ii[i+1] - ii[i]; 1650 aj = a->j + ii[i]; 1651 aa = a_a + ii[i]; 1652 sum = y[*ridx]; 1653 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1654 z[*ridx++] = sum; 1655 } 1656 } else { /* do not use compressed row format */ 1657 ii = a->i; 1658 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1659 aj = a->j; 1660 aa = a_a; 1661 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1662 #else 1663 for (i=0; i<m; i++) { 1664 n = ii[i+1] - ii[i]; 1665 aj = a->j + ii[i]; 1666 aa = a_a + ii[i]; 1667 sum = y[i]; 1668 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1669 z[i] = sum; 1670 } 1671 #endif 1672 } 1673 PetscCall(PetscLogFlops(2.0*a->nz)); 1674 PetscCall(VecRestoreArrayRead(xx,&x)); 1675 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1676 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1677 PetscFunctionReturn(0); 1678 } 1679 1680 /* 1681 Adds diagonal pointers to sparse matrix structure. 1682 */ 1683 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1684 { 1685 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1686 PetscInt i,j,m = A->rmap->n; 1687 PetscBool alreadySet = PETSC_TRUE; 1688 1689 PetscFunctionBegin; 1690 if (!a->diag) { 1691 PetscCall(PetscMalloc1(m,&a->diag)); 1692 PetscCall(PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt))); 1693 alreadySet = PETSC_FALSE; 1694 } 1695 for (i=0; i<A->rmap->n; i++) { 1696 /* If A's diagonal is already correctly set, this fast track enables cheap and repeated MatMarkDiagonal_SeqAIJ() calls */ 1697 if (alreadySet) { 1698 PetscInt pos = a->diag[i]; 1699 if (pos >= a->i[i] && pos < a->i[i+1] && a->j[pos] == i) continue; 1700 } 1701 1702 a->diag[i] = a->i[i+1]; 1703 for (j=a->i[i]; j<a->i[i+1]; j++) { 1704 if (a->j[j] == i) { 1705 a->diag[i] = j; 1706 break; 1707 } 1708 } 1709 } 1710 PetscFunctionReturn(0); 1711 } 1712 1713 PetscErrorCode MatShift_SeqAIJ(Mat A,PetscScalar v) 1714 { 1715 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1716 const PetscInt *diag = (const PetscInt*)a->diag; 1717 const PetscInt *ii = (const PetscInt*) a->i; 1718 PetscInt i,*mdiag = NULL; 1719 PetscInt cnt = 0; /* how many diagonals are missing */ 1720 1721 PetscFunctionBegin; 1722 if (!A->preallocated || !a->nz) { 1723 PetscCall(MatSeqAIJSetPreallocation(A,1,NULL)); 1724 PetscCall(MatShift_Basic(A,v)); 1725 PetscFunctionReturn(0); 1726 } 1727 1728 if (a->diagonaldense) { 1729 cnt = 0; 1730 } else { 1731 PetscCall(PetscCalloc1(A->rmap->n,&mdiag)); 1732 for (i=0; i<A->rmap->n; i++) { 1733 if (diag[i] >= ii[i+1]) { 1734 cnt++; 1735 mdiag[i] = 1; 1736 } 1737 } 1738 } 1739 if (!cnt) { 1740 PetscCall(MatShift_Basic(A,v)); 1741 } else { 1742 PetscScalar *olda = a->a; /* preserve pointers to current matrix nonzeros structure and values */ 1743 PetscInt *oldj = a->j, *oldi = a->i; 1744 PetscBool singlemalloc = a->singlemalloc,free_a = a->free_a,free_ij = a->free_ij; 1745 1746 a->a = NULL; 1747 a->j = NULL; 1748 a->i = NULL; 1749 /* increase the values in imax for each row where a diagonal is being inserted then reallocate the matrix data structures */ 1750 for (i=0; i<A->rmap->n; i++) { 1751 a->imax[i] += mdiag[i]; 1752 a->imax[i] = PetscMin(a->imax[i],A->cmap->n); 1753 } 1754 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax)); 1755 1756 /* copy old values into new matrix data structure */ 1757 for (i=0; i<A->rmap->n; i++) { 1758 PetscCall(MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES)); 1759 if (i < A->cmap->n) { 1760 PetscCall(MatSetValue(A,i,i,v,ADD_VALUES)); 1761 } 1762 } 1763 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1764 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1765 if (singlemalloc) { 1766 PetscCall(PetscFree3(olda,oldj,oldi)); 1767 } else { 1768 if (free_a) PetscCall(PetscFree(olda)); 1769 if (free_ij) PetscCall(PetscFree(oldj)); 1770 if (free_ij) PetscCall(PetscFree(oldi)); 1771 } 1772 } 1773 PetscCall(PetscFree(mdiag)); 1774 a->diagonaldense = PETSC_TRUE; 1775 PetscFunctionReturn(0); 1776 } 1777 1778 /* 1779 Checks for missing diagonals 1780 */ 1781 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1782 { 1783 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1784 PetscInt *diag,*ii = a->i,i; 1785 1786 PetscFunctionBegin; 1787 *missing = PETSC_FALSE; 1788 if (A->rmap->n > 0 && !ii) { 1789 *missing = PETSC_TRUE; 1790 if (d) *d = 0; 1791 PetscCall(PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n")); 1792 } else { 1793 PetscInt n; 1794 n = PetscMin(A->rmap->n, A->cmap->n); 1795 diag = a->diag; 1796 for (i=0; i<n; i++) { 1797 if (diag[i] >= ii[i+1]) { 1798 *missing = PETSC_TRUE; 1799 if (d) *d = i; 1800 PetscCall(PetscInfo(A,"Matrix is missing diagonal number %" PetscInt_FMT "\n",i)); 1801 break; 1802 } 1803 } 1804 } 1805 PetscFunctionReturn(0); 1806 } 1807 1808 #include <petscblaslapack.h> 1809 #include <petsc/private/kernels/blockinvert.h> 1810 1811 /* 1812 Note that values is allocated externally by the PC and then passed into this routine 1813 */ 1814 PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag) 1815 { 1816 PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots; 1817 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 1818 const PetscReal shift = 0.0; 1819 PetscInt ipvt[5]; 1820 PetscScalar work[25],*v_work; 1821 1822 PetscFunctionBegin; 1823 allowzeropivot = PetscNot(A->erroriffailure); 1824 for (i=0; i<nblocks; i++) ncnt += bsizes[i]; 1825 PetscCheck(ncnt == n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Total blocksizes %" PetscInt_FMT " doesn't match number matrix rows %" PetscInt_FMT,ncnt,n); 1826 for (i=0; i<nblocks; i++) { 1827 bsizemax = PetscMax(bsizemax,bsizes[i]); 1828 } 1829 PetscCall(PetscMalloc1(bsizemax,&indx)); 1830 if (bsizemax > 7) { 1831 PetscCall(PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots)); 1832 } 1833 ncnt = 0; 1834 for (i=0; i<nblocks; i++) { 1835 for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j; 1836 PetscCall(MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag)); 1837 switch (bsizes[i]) { 1838 case 1: 1839 *diag = 1.0/(*diag); 1840 break; 1841 case 2: 1842 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 1843 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1844 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 1845 break; 1846 case 3: 1847 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 1848 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1849 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 1850 break; 1851 case 4: 1852 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 1853 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1854 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 1855 break; 1856 case 5: 1857 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 1858 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1859 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 1860 break; 1861 case 6: 1862 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 1863 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1864 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 1865 break; 1866 case 7: 1867 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 1868 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1869 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 1870 break; 1871 default: 1872 PetscCall(PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 1873 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1874 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bsizes[i])); 1875 } 1876 ncnt += bsizes[i]; 1877 diag += bsizes[i]*bsizes[i]; 1878 } 1879 if (bsizemax > 7) { 1880 PetscCall(PetscFree2(v_work,v_pivots)); 1881 } 1882 PetscCall(PetscFree(indx)); 1883 PetscFunctionReturn(0); 1884 } 1885 1886 /* 1887 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1888 */ 1889 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1890 { 1891 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1892 PetscInt i,*diag,m = A->rmap->n; 1893 const MatScalar *v; 1894 PetscScalar *idiag,*mdiag; 1895 1896 PetscFunctionBegin; 1897 if (a->idiagvalid) PetscFunctionReturn(0); 1898 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1899 diag = a->diag; 1900 if (!a->idiag) { 1901 PetscCall(PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work)); 1902 PetscCall(PetscLogObjectMemory((PetscObject)A,3*m*sizeof(PetscScalar))); 1903 } 1904 1905 mdiag = a->mdiag; 1906 idiag = a->idiag; 1907 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 1908 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1909 for (i=0; i<m; i++) { 1910 mdiag[i] = v[diag[i]]; 1911 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1912 if (PetscRealPart(fshift)) { 1913 PetscCall(PetscInfo(A,"Zero diagonal on row %" PetscInt_FMT "\n",i)); 1914 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1915 A->factorerror_zeropivot_value = 0.0; 1916 A->factorerror_zeropivot_row = i; 1917 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %" PetscInt_FMT,i); 1918 } 1919 idiag[i] = 1.0/v[diag[i]]; 1920 } 1921 PetscCall(PetscLogFlops(m)); 1922 } else { 1923 for (i=0; i<m; i++) { 1924 mdiag[i] = v[diag[i]]; 1925 idiag[i] = omega/(fshift + v[diag[i]]); 1926 } 1927 PetscCall(PetscLogFlops(2.0*m)); 1928 } 1929 a->idiagvalid = PETSC_TRUE; 1930 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 1931 PetscFunctionReturn(0); 1932 } 1933 1934 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1935 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1936 { 1937 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1938 PetscScalar *x,d,sum,*t,scale; 1939 const MatScalar *v,*idiag=NULL,*mdiag,*aa; 1940 const PetscScalar *b, *bs,*xb, *ts; 1941 PetscInt n,m = A->rmap->n,i; 1942 const PetscInt *idx,*diag; 1943 1944 PetscFunctionBegin; 1945 if (a->inode.use && a->inode.checked && omega == 1.0 && fshift == 0.0) { 1946 PetscCall(MatSOR_SeqAIJ_Inode(A,bb,omega,flag,fshift,its,lits,xx)); 1947 PetscFunctionReturn(0); 1948 } 1949 its = its*lits; 1950 1951 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1952 if (!a->idiagvalid) PetscCall(MatInvertDiagonal_SeqAIJ(A,omega,fshift)); 1953 a->fshift = fshift; 1954 a->omega = omega; 1955 1956 diag = a->diag; 1957 t = a->ssor_work; 1958 idiag = a->idiag; 1959 mdiag = a->mdiag; 1960 1961 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1962 PetscCall(VecGetArray(xx,&x)); 1963 PetscCall(VecGetArrayRead(bb,&b)); 1964 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1965 if (flag == SOR_APPLY_UPPER) { 1966 /* apply (U + D/omega) to the vector */ 1967 bs = b; 1968 for (i=0; i<m; i++) { 1969 d = fshift + mdiag[i]; 1970 n = a->i[i+1] - diag[i] - 1; 1971 idx = a->j + diag[i] + 1; 1972 v = aa + diag[i] + 1; 1973 sum = b[i]*d/omega; 1974 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1975 x[i] = sum; 1976 } 1977 PetscCall(VecRestoreArray(xx,&x)); 1978 PetscCall(VecRestoreArrayRead(bb,&b)); 1979 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1980 PetscCall(PetscLogFlops(a->nz)); 1981 PetscFunctionReturn(0); 1982 } 1983 1984 PetscCheck(flag != SOR_APPLY_LOWER,PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1985 else if (flag & SOR_EISENSTAT) { 1986 /* Let A = L + U + D; where L is lower triangular, 1987 U is upper triangular, E = D/omega; This routine applies 1988 1989 (L + E)^{-1} A (U + E)^{-1} 1990 1991 to a vector efficiently using Eisenstat's trick. 1992 */ 1993 scale = (2.0/omega) - 1.0; 1994 1995 /* x = (E + U)^{-1} b */ 1996 for (i=m-1; i>=0; i--) { 1997 n = a->i[i+1] - diag[i] - 1; 1998 idx = a->j + diag[i] + 1; 1999 v = aa + diag[i] + 1; 2000 sum = b[i]; 2001 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2002 x[i] = sum*idiag[i]; 2003 } 2004 2005 /* t = b - (2*E - D)x */ 2006 v = aa; 2007 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 2008 2009 /* t = (E + L)^{-1}t */ 2010 ts = t; 2011 diag = a->diag; 2012 for (i=0; i<m; i++) { 2013 n = diag[i] - a->i[i]; 2014 idx = a->j + a->i[i]; 2015 v = aa + a->i[i]; 2016 sum = t[i]; 2017 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 2018 t[i] = sum*idiag[i]; 2019 /* x = x + t */ 2020 x[i] += t[i]; 2021 } 2022 2023 PetscCall(PetscLogFlops(6.0*m-1 + 2.0*a->nz)); 2024 PetscCall(VecRestoreArray(xx,&x)); 2025 PetscCall(VecRestoreArrayRead(bb,&b)); 2026 PetscFunctionReturn(0); 2027 } 2028 if (flag & SOR_ZERO_INITIAL_GUESS) { 2029 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2030 for (i=0; i<m; i++) { 2031 n = diag[i] - a->i[i]; 2032 idx = a->j + a->i[i]; 2033 v = aa + a->i[i]; 2034 sum = b[i]; 2035 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2036 t[i] = sum; 2037 x[i] = sum*idiag[i]; 2038 } 2039 xb = t; 2040 PetscCall(PetscLogFlops(a->nz)); 2041 } else xb = b; 2042 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2043 for (i=m-1; i>=0; i--) { 2044 n = a->i[i+1] - diag[i] - 1; 2045 idx = a->j + diag[i] + 1; 2046 v = aa + diag[i] + 1; 2047 sum = xb[i]; 2048 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2049 if (xb == b) { 2050 x[i] = sum*idiag[i]; 2051 } else { 2052 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2053 } 2054 } 2055 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2056 } 2057 its--; 2058 } 2059 while (its--) { 2060 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2061 for (i=0; i<m; i++) { 2062 /* lower */ 2063 n = diag[i] - a->i[i]; 2064 idx = a->j + a->i[i]; 2065 v = aa + a->i[i]; 2066 sum = b[i]; 2067 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2068 t[i] = sum; /* save application of the lower-triangular part */ 2069 /* upper */ 2070 n = a->i[i+1] - diag[i] - 1; 2071 idx = a->j + diag[i] + 1; 2072 v = aa + diag[i] + 1; 2073 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2074 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2075 } 2076 xb = t; 2077 PetscCall(PetscLogFlops(2.0*a->nz)); 2078 } else xb = b; 2079 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2080 for (i=m-1; i>=0; i--) { 2081 sum = xb[i]; 2082 if (xb == b) { 2083 /* whole matrix (no checkpointing available) */ 2084 n = a->i[i+1] - a->i[i]; 2085 idx = a->j + a->i[i]; 2086 v = aa + a->i[i]; 2087 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2088 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 2089 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 2090 n = a->i[i+1] - diag[i] - 1; 2091 idx = a->j + diag[i] + 1; 2092 v = aa + diag[i] + 1; 2093 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2094 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2095 } 2096 } 2097 if (xb == b) { 2098 PetscCall(PetscLogFlops(2.0*a->nz)); 2099 } else { 2100 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2101 } 2102 } 2103 } 2104 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2105 PetscCall(VecRestoreArray(xx,&x)); 2106 PetscCall(VecRestoreArrayRead(bb,&b)); 2107 PetscFunctionReturn(0); 2108 } 2109 2110 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 2111 { 2112 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2113 2114 PetscFunctionBegin; 2115 info->block_size = 1.0; 2116 info->nz_allocated = a->maxnz; 2117 info->nz_used = a->nz; 2118 info->nz_unneeded = (a->maxnz - a->nz); 2119 info->assemblies = A->num_ass; 2120 info->mallocs = A->info.mallocs; 2121 info->memory = ((PetscObject)A)->mem; 2122 if (A->factortype) { 2123 info->fill_ratio_given = A->info.fill_ratio_given; 2124 info->fill_ratio_needed = A->info.fill_ratio_needed; 2125 info->factor_mallocs = A->info.factor_mallocs; 2126 } else { 2127 info->fill_ratio_given = 0; 2128 info->fill_ratio_needed = 0; 2129 info->factor_mallocs = 0; 2130 } 2131 PetscFunctionReturn(0); 2132 } 2133 2134 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2135 { 2136 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2137 PetscInt i,m = A->rmap->n - 1; 2138 const PetscScalar *xx; 2139 PetscScalar *bb,*aa; 2140 PetscInt d = 0; 2141 2142 PetscFunctionBegin; 2143 if (x && b) { 2144 PetscCall(VecGetArrayRead(x,&xx)); 2145 PetscCall(VecGetArray(b,&bb)); 2146 for (i=0; i<N; i++) { 2147 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2148 if (rows[i] >= A->cmap->n) continue; 2149 bb[rows[i]] = diag*xx[rows[i]]; 2150 } 2151 PetscCall(VecRestoreArrayRead(x,&xx)); 2152 PetscCall(VecRestoreArray(b,&bb)); 2153 } 2154 2155 PetscCall(MatSeqAIJGetArray(A,&aa)); 2156 if (a->keepnonzeropattern) { 2157 for (i=0; i<N; i++) { 2158 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2159 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2160 } 2161 if (diag != 0.0) { 2162 for (i=0; i<N; i++) { 2163 d = rows[i]; 2164 if (rows[i] >= A->cmap->n) continue; 2165 PetscCheck(a->diag[d] < a->i[d+1],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in the zeroed row %" PetscInt_FMT,d); 2166 } 2167 for (i=0; i<N; i++) { 2168 if (rows[i] >= A->cmap->n) continue; 2169 aa[a->diag[rows[i]]] = diag; 2170 } 2171 } 2172 } else { 2173 if (diag != 0.0) { 2174 for (i=0; i<N; i++) { 2175 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2176 if (a->ilen[rows[i]] > 0) { 2177 if (rows[i] >= A->cmap->n) { 2178 a->ilen[rows[i]] = 0; 2179 } else { 2180 a->ilen[rows[i]] = 1; 2181 aa[a->i[rows[i]]] = diag; 2182 a->j[a->i[rows[i]]] = rows[i]; 2183 } 2184 } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */ 2185 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2186 } 2187 } 2188 } else { 2189 for (i=0; i<N; i++) { 2190 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2191 a->ilen[rows[i]] = 0; 2192 } 2193 } 2194 A->nonzerostate++; 2195 } 2196 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2197 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2198 PetscFunctionReturn(0); 2199 } 2200 2201 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2202 { 2203 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2204 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2205 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2206 const PetscScalar *xx; 2207 PetscScalar *bb,*aa; 2208 2209 PetscFunctionBegin; 2210 if (!N) PetscFunctionReturn(0); 2211 PetscCall(MatSeqAIJGetArray(A,&aa)); 2212 if (x && b) { 2213 PetscCall(VecGetArrayRead(x,&xx)); 2214 PetscCall(VecGetArray(b,&bb)); 2215 vecs = PETSC_TRUE; 2216 } 2217 PetscCall(PetscCalloc1(A->rmap->n,&zeroed)); 2218 for (i=0; i<N; i++) { 2219 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2220 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2221 2222 zeroed[rows[i]] = PETSC_TRUE; 2223 } 2224 for (i=0; i<A->rmap->n; i++) { 2225 if (!zeroed[i]) { 2226 for (j=a->i[i]; j<a->i[i+1]; j++) { 2227 if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) { 2228 if (vecs) bb[i] -= aa[j]*xx[a->j[j]]; 2229 aa[j] = 0.0; 2230 } 2231 } 2232 } else if (vecs && i < A->cmap->N) bb[i] = diag*xx[i]; 2233 } 2234 if (x && b) { 2235 PetscCall(VecRestoreArrayRead(x,&xx)); 2236 PetscCall(VecRestoreArray(b,&bb)); 2237 } 2238 PetscCall(PetscFree(zeroed)); 2239 if (diag != 0.0) { 2240 PetscCall(MatMissingDiagonal_SeqAIJ(A,&missing,&d)); 2241 if (missing) { 2242 for (i=0; i<N; i++) { 2243 if (rows[i] >= A->cmap->N) continue; 2244 PetscCheckFalse(a->nonew && rows[i] >= d,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %" PetscInt_FMT " (%" PetscInt_FMT ")",d,rows[i]); 2245 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2246 } 2247 } else { 2248 for (i=0; i<N; i++) { 2249 aa[a->diag[rows[i]]] = diag; 2250 } 2251 } 2252 } 2253 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2254 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2255 PetscFunctionReturn(0); 2256 } 2257 2258 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2259 { 2260 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2261 const PetscScalar *aa; 2262 PetscInt *itmp; 2263 2264 PetscFunctionBegin; 2265 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2266 *nz = a->i[row+1] - a->i[row]; 2267 if (v) *v = (PetscScalar*)(aa + a->i[row]); 2268 if (idx) { 2269 itmp = a->j + a->i[row]; 2270 if (*nz) *idx = itmp; 2271 else *idx = NULL; 2272 } 2273 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2274 PetscFunctionReturn(0); 2275 } 2276 2277 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2278 { 2279 PetscFunctionBegin; 2280 if (nz) *nz = 0; 2281 if (idx) *idx = NULL; 2282 if (v) *v = NULL; 2283 PetscFunctionReturn(0); 2284 } 2285 2286 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2287 { 2288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2289 const MatScalar *v; 2290 PetscReal sum = 0.0; 2291 PetscInt i,j; 2292 2293 PetscFunctionBegin; 2294 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 2295 if (type == NORM_FROBENIUS) { 2296 #if defined(PETSC_USE_REAL___FP16) 2297 PetscBLASInt one = 1,nz = a->nz; 2298 PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&nz,v,&one)); 2299 #else 2300 for (i=0; i<a->nz; i++) { 2301 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2302 } 2303 *nrm = PetscSqrtReal(sum); 2304 #endif 2305 PetscCall(PetscLogFlops(2.0*a->nz)); 2306 } else if (type == NORM_1) { 2307 PetscReal *tmp; 2308 PetscInt *jj = a->j; 2309 PetscCall(PetscCalloc1(A->cmap->n+1,&tmp)); 2310 *nrm = 0.0; 2311 for (j=0; j<a->nz; j++) { 2312 tmp[*jj++] += PetscAbsScalar(*v); v++; 2313 } 2314 for (j=0; j<A->cmap->n; j++) { 2315 if (tmp[j] > *nrm) *nrm = tmp[j]; 2316 } 2317 PetscCall(PetscFree(tmp)); 2318 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2319 } else if (type == NORM_INFINITY) { 2320 *nrm = 0.0; 2321 for (j=0; j<A->rmap->n; j++) { 2322 const PetscScalar *v2 = v + a->i[j]; 2323 sum = 0.0; 2324 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2325 sum += PetscAbsScalar(*v2); v2++; 2326 } 2327 if (sum > *nrm) *nrm = sum; 2328 } 2329 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2330 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2331 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 2332 PetscFunctionReturn(0); 2333 } 2334 2335 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2336 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2337 { 2338 PetscInt i,j,anzj; 2339 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2340 PetscInt an=A->cmap->N,am=A->rmap->N; 2341 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2342 2343 PetscFunctionBegin; 2344 /* Allocate space for symbolic transpose info and work array */ 2345 PetscCall(PetscCalloc1(an+1,&ati)); 2346 PetscCall(PetscMalloc1(ai[am],&atj)); 2347 PetscCall(PetscMalloc1(an,&atfill)); 2348 2349 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2350 /* Note: offset by 1 for fast conversion into csr format. */ 2351 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2352 /* Form ati for csr format of A^T. */ 2353 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2354 2355 /* Copy ati into atfill so we have locations of the next free space in atj */ 2356 PetscCall(PetscArraycpy(atfill,ati,an)); 2357 2358 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2359 for (i=0;i<am;i++) { 2360 anzj = ai[i+1] - ai[i]; 2361 for (j=0;j<anzj;j++) { 2362 atj[atfill[*aj]] = i; 2363 atfill[*aj++] += 1; 2364 } 2365 } 2366 2367 /* Clean up temporary space and complete requests. */ 2368 PetscCall(PetscFree(atfill)); 2369 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B)); 2370 PetscCall(MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs))); 2371 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 2372 2373 b = (Mat_SeqAIJ*)((*B)->data); 2374 b->free_a = PETSC_FALSE; 2375 b->free_ij = PETSC_TRUE; 2376 b->nonew = 0; 2377 PetscFunctionReturn(0); 2378 } 2379 2380 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2381 { 2382 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2383 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2384 const MatScalar *va,*vb; 2385 PetscInt ma,na,mb,nb, i; 2386 2387 PetscFunctionBegin; 2388 PetscCall(MatGetSize(A,&ma,&na)); 2389 PetscCall(MatGetSize(B,&mb,&nb)); 2390 if (ma!=nb || na!=mb) { 2391 *f = PETSC_FALSE; 2392 PetscFunctionReturn(0); 2393 } 2394 PetscCall(MatSeqAIJGetArrayRead(A,&va)); 2395 PetscCall(MatSeqAIJGetArrayRead(B,&vb)); 2396 aii = aij->i; bii = bij->i; 2397 adx = aij->j; bdx = bij->j; 2398 PetscCall(PetscMalloc1(ma,&aptr)); 2399 PetscCall(PetscMalloc1(mb,&bptr)); 2400 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2401 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2402 2403 *f = PETSC_TRUE; 2404 for (i=0; i<ma; i++) { 2405 while (aptr[i]<aii[i+1]) { 2406 PetscInt idc,idr; 2407 PetscScalar vc,vr; 2408 /* column/row index/value */ 2409 idc = adx[aptr[i]]; 2410 idr = bdx[bptr[idc]]; 2411 vc = va[aptr[i]]; 2412 vr = vb[bptr[idc]]; 2413 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2414 *f = PETSC_FALSE; 2415 goto done; 2416 } else { 2417 aptr[i]++; 2418 if (B || i!=idc) bptr[idc]++; 2419 } 2420 } 2421 } 2422 done: 2423 PetscCall(PetscFree(aptr)); 2424 PetscCall(PetscFree(bptr)); 2425 PetscCall(MatSeqAIJRestoreArrayRead(A,&va)); 2426 PetscCall(MatSeqAIJRestoreArrayRead(B,&vb)); 2427 PetscFunctionReturn(0); 2428 } 2429 2430 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2431 { 2432 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2433 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2434 MatScalar *va,*vb; 2435 PetscInt ma,na,mb,nb, i; 2436 2437 PetscFunctionBegin; 2438 PetscCall(MatGetSize(A,&ma,&na)); 2439 PetscCall(MatGetSize(B,&mb,&nb)); 2440 if (ma!=nb || na!=mb) { 2441 *f = PETSC_FALSE; 2442 PetscFunctionReturn(0); 2443 } 2444 aii = aij->i; bii = bij->i; 2445 adx = aij->j; bdx = bij->j; 2446 va = aij->a; vb = bij->a; 2447 PetscCall(PetscMalloc1(ma,&aptr)); 2448 PetscCall(PetscMalloc1(mb,&bptr)); 2449 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2450 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2451 2452 *f = PETSC_TRUE; 2453 for (i=0; i<ma; i++) { 2454 while (aptr[i]<aii[i+1]) { 2455 PetscInt idc,idr; 2456 PetscScalar vc,vr; 2457 /* column/row index/value */ 2458 idc = adx[aptr[i]]; 2459 idr = bdx[bptr[idc]]; 2460 vc = va[aptr[i]]; 2461 vr = vb[bptr[idc]]; 2462 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2463 *f = PETSC_FALSE; 2464 goto done; 2465 } else { 2466 aptr[i]++; 2467 if (B || i!=idc) bptr[idc]++; 2468 } 2469 } 2470 } 2471 done: 2472 PetscCall(PetscFree(aptr)); 2473 PetscCall(PetscFree(bptr)); 2474 PetscFunctionReturn(0); 2475 } 2476 2477 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2478 { 2479 PetscFunctionBegin; 2480 PetscCall(MatIsTranspose_SeqAIJ(A,A,tol,f)); 2481 PetscFunctionReturn(0); 2482 } 2483 2484 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2485 { 2486 PetscFunctionBegin; 2487 PetscCall(MatIsHermitianTranspose_SeqAIJ(A,A,tol,f)); 2488 PetscFunctionReturn(0); 2489 } 2490 2491 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2492 { 2493 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2494 const PetscScalar *l,*r; 2495 PetscScalar x; 2496 MatScalar *v; 2497 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2498 const PetscInt *jj; 2499 2500 PetscFunctionBegin; 2501 if (ll) { 2502 /* The local size is used so that VecMPI can be passed to this routine 2503 by MatDiagonalScale_MPIAIJ */ 2504 PetscCall(VecGetLocalSize(ll,&m)); 2505 PetscCheck(m == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2506 PetscCall(VecGetArrayRead(ll,&l)); 2507 PetscCall(MatSeqAIJGetArray(A,&v)); 2508 for (i=0; i<m; i++) { 2509 x = l[i]; 2510 M = a->i[i+1] - a->i[i]; 2511 for (j=0; j<M; j++) (*v++) *= x; 2512 } 2513 PetscCall(VecRestoreArrayRead(ll,&l)); 2514 PetscCall(PetscLogFlops(nz)); 2515 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2516 } 2517 if (rr) { 2518 PetscCall(VecGetLocalSize(rr,&n)); 2519 PetscCheck(n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2520 PetscCall(VecGetArrayRead(rr,&r)); 2521 PetscCall(MatSeqAIJGetArray(A,&v)); 2522 jj = a->j; 2523 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2524 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2525 PetscCall(VecRestoreArrayRead(rr,&r)); 2526 PetscCall(PetscLogFlops(nz)); 2527 } 2528 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 2529 PetscFunctionReturn(0); 2530 } 2531 2532 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2533 { 2534 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2535 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2536 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2537 const PetscInt *irow,*icol; 2538 const PetscScalar *aa; 2539 PetscInt nrows,ncols; 2540 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2541 MatScalar *a_new,*mat_a; 2542 Mat C; 2543 PetscBool stride; 2544 2545 PetscFunctionBegin; 2546 PetscCall(ISGetIndices(isrow,&irow)); 2547 PetscCall(ISGetLocalSize(isrow,&nrows)); 2548 PetscCall(ISGetLocalSize(iscol,&ncols)); 2549 2550 PetscCall(PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride)); 2551 if (stride) { 2552 PetscCall(ISStrideGetInfo(iscol,&first,&step)); 2553 } else { 2554 first = 0; 2555 step = 0; 2556 } 2557 if (stride && step == 1) { 2558 /* special case of contiguous rows */ 2559 PetscCall(PetscMalloc2(nrows,&lens,nrows,&starts)); 2560 /* loop over new rows determining lens and starting points */ 2561 for (i=0; i<nrows; i++) { 2562 kstart = ai[irow[i]]; 2563 kend = kstart + ailen[irow[i]]; 2564 starts[i] = kstart; 2565 for (k=kstart; k<kend; k++) { 2566 if (aj[k] >= first) { 2567 starts[i] = k; 2568 break; 2569 } 2570 } 2571 sum = 0; 2572 while (k < kend) { 2573 if (aj[k++] >= first+ncols) break; 2574 sum++; 2575 } 2576 lens[i] = sum; 2577 } 2578 /* create submatrix */ 2579 if (scall == MAT_REUSE_MATRIX) { 2580 PetscInt n_cols,n_rows; 2581 PetscCall(MatGetSize(*B,&n_rows,&n_cols)); 2582 PetscCheckFalse(n_rows != nrows || n_cols != ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2583 PetscCall(MatZeroEntries(*B)); 2584 C = *B; 2585 } else { 2586 PetscInt rbs,cbs; 2587 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2588 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2589 PetscCall(ISGetBlockSize(isrow,&rbs)); 2590 PetscCall(ISGetBlockSize(iscol,&cbs)); 2591 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2592 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2593 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2594 } 2595 c = (Mat_SeqAIJ*)C->data; 2596 2597 /* loop over rows inserting into submatrix */ 2598 a_new = c->a; 2599 j_new = c->j; 2600 i_new = c->i; 2601 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2602 for (i=0; i<nrows; i++) { 2603 ii = starts[i]; 2604 lensi = lens[i]; 2605 for (k=0; k<lensi; k++) { 2606 *j_new++ = aj[ii+k] - first; 2607 } 2608 PetscCall(PetscArraycpy(a_new,aa + starts[i],lensi)); 2609 a_new += lensi; 2610 i_new[i+1] = i_new[i] + lensi; 2611 c->ilen[i] = lensi; 2612 } 2613 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2614 PetscCall(PetscFree2(lens,starts)); 2615 } else { 2616 PetscCall(ISGetIndices(iscol,&icol)); 2617 PetscCall(PetscCalloc1(oldcols,&smap)); 2618 PetscCall(PetscMalloc1(1+nrows,&lens)); 2619 for (i=0; i<ncols; i++) { 2620 PetscCheck(icol[i] < oldcols,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requesting column beyond largest column icol[%" PetscInt_FMT "] %" PetscInt_FMT " >= A->cmap->n %" PetscInt_FMT,i,icol[i],oldcols); 2621 smap[icol[i]] = i+1; 2622 } 2623 2624 /* determine lens of each row */ 2625 for (i=0; i<nrows; i++) { 2626 kstart = ai[irow[i]]; 2627 kend = kstart + a->ilen[irow[i]]; 2628 lens[i] = 0; 2629 for (k=kstart; k<kend; k++) { 2630 if (smap[aj[k]]) { 2631 lens[i]++; 2632 } 2633 } 2634 } 2635 /* Create and fill new matrix */ 2636 if (scall == MAT_REUSE_MATRIX) { 2637 PetscBool equal; 2638 2639 c = (Mat_SeqAIJ*)((*B)->data); 2640 PetscCheckFalse((*B)->rmap->n != nrows || (*B)->cmap->n != ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2641 PetscCall(PetscArraycmp(c->ilen,lens,(*B)->rmap->n,&equal)); 2642 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2643 PetscCall(PetscArrayzero(c->ilen,(*B)->rmap->n)); 2644 C = *B; 2645 } else { 2646 PetscInt rbs,cbs; 2647 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2648 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2649 PetscCall(ISGetBlockSize(isrow,&rbs)); 2650 PetscCall(ISGetBlockSize(iscol,&cbs)); 2651 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2652 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2653 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2654 } 2655 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2656 c = (Mat_SeqAIJ*)(C->data); 2657 for (i=0; i<nrows; i++) { 2658 row = irow[i]; 2659 kstart = ai[row]; 2660 kend = kstart + a->ilen[row]; 2661 mat_i = c->i[i]; 2662 mat_j = c->j + mat_i; 2663 mat_a = c->a + mat_i; 2664 mat_ilen = c->ilen + i; 2665 for (k=kstart; k<kend; k++) { 2666 if ((tcol=smap[a->j[k]])) { 2667 *mat_j++ = tcol - 1; 2668 *mat_a++ = aa[k]; 2669 (*mat_ilen)++; 2670 2671 } 2672 } 2673 } 2674 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2675 /* Free work space */ 2676 PetscCall(ISRestoreIndices(iscol,&icol)); 2677 PetscCall(PetscFree(smap)); 2678 PetscCall(PetscFree(lens)); 2679 /* sort */ 2680 for (i = 0; i < nrows; i++) { 2681 PetscInt ilen; 2682 2683 mat_i = c->i[i]; 2684 mat_j = c->j + mat_i; 2685 mat_a = c->a + mat_i; 2686 ilen = c->ilen[i]; 2687 PetscCall(PetscSortIntWithScalarArray(ilen,mat_j,mat_a)); 2688 } 2689 } 2690 #if defined(PETSC_HAVE_DEVICE) 2691 PetscCall(MatBindToCPU(C,A->boundtocpu)); 2692 #endif 2693 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 2694 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 2695 2696 PetscCall(ISRestoreIndices(isrow,&irow)); 2697 *B = C; 2698 PetscFunctionReturn(0); 2699 } 2700 2701 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2702 { 2703 Mat B; 2704 2705 PetscFunctionBegin; 2706 if (scall == MAT_INITIAL_MATRIX) { 2707 PetscCall(MatCreate(subComm,&B)); 2708 PetscCall(MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n)); 2709 PetscCall(MatSetBlockSizesFromMats(B,mat,mat)); 2710 PetscCall(MatSetType(B,MATSEQAIJ)); 2711 PetscCall(MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE)); 2712 *subMat = B; 2713 } else { 2714 PetscCall(MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN)); 2715 } 2716 PetscFunctionReturn(0); 2717 } 2718 2719 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2720 { 2721 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2722 Mat outA; 2723 PetscBool row_identity,col_identity; 2724 2725 PetscFunctionBegin; 2726 PetscCheck(info->levels == 0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2727 2728 PetscCall(ISIdentity(row,&row_identity)); 2729 PetscCall(ISIdentity(col,&col_identity)); 2730 2731 outA = inA; 2732 outA->factortype = MAT_FACTOR_LU; 2733 PetscCall(PetscFree(inA->solvertype)); 2734 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype)); 2735 2736 PetscCall(PetscObjectReference((PetscObject)row)); 2737 PetscCall(ISDestroy(&a->row)); 2738 2739 a->row = row; 2740 2741 PetscCall(PetscObjectReference((PetscObject)col)); 2742 PetscCall(ISDestroy(&a->col)); 2743 2744 a->col = col; 2745 2746 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2747 PetscCall(ISDestroy(&a->icol)); 2748 PetscCall(ISInvertPermutation(col,PETSC_DECIDE,&a->icol)); 2749 PetscCall(PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol)); 2750 2751 if (!a->solve_work) { /* this matrix may have been factored before */ 2752 PetscCall(PetscMalloc1(inA->rmap->n+1,&a->solve_work)); 2753 PetscCall(PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar))); 2754 } 2755 2756 PetscCall(MatMarkDiagonal_SeqAIJ(inA)); 2757 if (row_identity && col_identity) { 2758 PetscCall(MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info)); 2759 } else { 2760 PetscCall(MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info)); 2761 } 2762 PetscFunctionReturn(0); 2763 } 2764 2765 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2766 { 2767 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2768 PetscScalar *v; 2769 PetscBLASInt one = 1,bnz; 2770 2771 PetscFunctionBegin; 2772 PetscCall(MatSeqAIJGetArray(inA,&v)); 2773 PetscCall(PetscBLASIntCast(a->nz,&bnz)); 2774 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&alpha,v,&one)); 2775 PetscCall(PetscLogFlops(a->nz)); 2776 PetscCall(MatSeqAIJRestoreArray(inA,&v)); 2777 PetscCall(MatSeqAIJInvalidateDiagonal(inA)); 2778 PetscFunctionReturn(0); 2779 } 2780 2781 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2782 { 2783 PetscInt i; 2784 2785 PetscFunctionBegin; 2786 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2787 PetscCall(PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr)); 2788 2789 for (i=0; i<submatj->nrqr; ++i) { 2790 PetscCall(PetscFree(submatj->sbuf2[i])); 2791 } 2792 PetscCall(PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1)); 2793 2794 if (submatj->rbuf1) { 2795 PetscCall(PetscFree(submatj->rbuf1[0])); 2796 PetscCall(PetscFree(submatj->rbuf1)); 2797 } 2798 2799 for (i=0; i<submatj->nrqs; ++i) { 2800 PetscCall(PetscFree(submatj->rbuf3[i])); 2801 } 2802 PetscCall(PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3)); 2803 PetscCall(PetscFree(submatj->pa)); 2804 } 2805 2806 #if defined(PETSC_USE_CTABLE) 2807 PetscCall(PetscTableDestroy((PetscTable*)&submatj->rmap)); 2808 if (submatj->cmap_loc) PetscCall(PetscFree(submatj->cmap_loc)); 2809 PetscCall(PetscFree(submatj->rmap_loc)); 2810 #else 2811 PetscCall(PetscFree(submatj->rmap)); 2812 #endif 2813 2814 if (!submatj->allcolumns) { 2815 #if defined(PETSC_USE_CTABLE) 2816 PetscCall(PetscTableDestroy((PetscTable*)&submatj->cmap)); 2817 #else 2818 PetscCall(PetscFree(submatj->cmap)); 2819 #endif 2820 } 2821 PetscCall(PetscFree(submatj->row2proc)); 2822 2823 PetscCall(PetscFree(submatj)); 2824 PetscFunctionReturn(0); 2825 } 2826 2827 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2828 { 2829 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2830 Mat_SubSppt *submatj = c->submatis1; 2831 2832 PetscFunctionBegin; 2833 PetscCall((*submatj->destroy)(C)); 2834 PetscCall(MatDestroySubMatrix_Private(submatj)); 2835 PetscFunctionReturn(0); 2836 } 2837 2838 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2839 { 2840 PetscInt i; 2841 Mat C; 2842 Mat_SeqAIJ *c; 2843 Mat_SubSppt *submatj; 2844 2845 PetscFunctionBegin; 2846 for (i=0; i<n; i++) { 2847 C = (*mat)[i]; 2848 c = (Mat_SeqAIJ*)C->data; 2849 submatj = c->submatis1; 2850 if (submatj) { 2851 if (--((PetscObject)C)->refct <= 0) { 2852 PetscCall((*submatj->destroy)(C)); 2853 PetscCall(MatDestroySubMatrix_Private(submatj)); 2854 PetscCall(PetscFree(C->defaultvectype)); 2855 PetscCall(PetscLayoutDestroy(&C->rmap)); 2856 PetscCall(PetscLayoutDestroy(&C->cmap)); 2857 PetscCall(PetscHeaderDestroy(&C)); 2858 } 2859 } else { 2860 PetscCall(MatDestroy(&C)); 2861 } 2862 } 2863 2864 /* Destroy Dummy submatrices created for reuse */ 2865 PetscCall(MatDestroySubMatrices_Dummy(n,mat)); 2866 2867 PetscCall(PetscFree(*mat)); 2868 PetscFunctionReturn(0); 2869 } 2870 2871 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2872 { 2873 PetscInt i; 2874 2875 PetscFunctionBegin; 2876 if (scall == MAT_INITIAL_MATRIX) { 2877 PetscCall(PetscCalloc1(n+1,B)); 2878 } 2879 2880 for (i=0; i<n; i++) { 2881 PetscCall(MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i])); 2882 } 2883 PetscFunctionReturn(0); 2884 } 2885 2886 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2887 { 2888 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2889 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2890 const PetscInt *idx; 2891 PetscInt start,end,*ai,*aj; 2892 PetscBT table; 2893 2894 PetscFunctionBegin; 2895 m = A->rmap->n; 2896 ai = a->i; 2897 aj = a->j; 2898 2899 PetscCheck(ov >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2900 2901 PetscCall(PetscMalloc1(m+1,&nidx)); 2902 PetscCall(PetscBTCreate(m,&table)); 2903 2904 for (i=0; i<is_max; i++) { 2905 /* Initialize the two local arrays */ 2906 isz = 0; 2907 PetscCall(PetscBTMemzero(m,table)); 2908 2909 /* Extract the indices, assume there can be duplicate entries */ 2910 PetscCall(ISGetIndices(is[i],&idx)); 2911 PetscCall(ISGetLocalSize(is[i],&n)); 2912 2913 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2914 for (j=0; j<n; ++j) { 2915 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2916 } 2917 PetscCall(ISRestoreIndices(is[i],&idx)); 2918 PetscCall(ISDestroy(&is[i])); 2919 2920 k = 0; 2921 for (j=0; j<ov; j++) { /* for each overlap */ 2922 n = isz; 2923 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2924 row = nidx[k]; 2925 start = ai[row]; 2926 end = ai[row+1]; 2927 for (l = start; l<end; l++) { 2928 val = aj[l]; 2929 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2930 } 2931 } 2932 } 2933 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i))); 2934 } 2935 PetscCall(PetscBTDestroy(&table)); 2936 PetscCall(PetscFree(nidx)); 2937 PetscFunctionReturn(0); 2938 } 2939 2940 /* -------------------------------------------------------------- */ 2941 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2942 { 2943 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2944 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2945 const PetscInt *row,*col; 2946 PetscInt *cnew,j,*lens; 2947 IS icolp,irowp; 2948 PetscInt *cwork = NULL; 2949 PetscScalar *vwork = NULL; 2950 2951 PetscFunctionBegin; 2952 PetscCall(ISInvertPermutation(rowp,PETSC_DECIDE,&irowp)); 2953 PetscCall(ISGetIndices(irowp,&row)); 2954 PetscCall(ISInvertPermutation(colp,PETSC_DECIDE,&icolp)); 2955 PetscCall(ISGetIndices(icolp,&col)); 2956 2957 /* determine lengths of permuted rows */ 2958 PetscCall(PetscMalloc1(m+1,&lens)); 2959 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2960 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 2961 PetscCall(MatSetSizes(*B,m,n,m,n)); 2962 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 2963 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 2964 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens)); 2965 PetscCall(PetscFree(lens)); 2966 2967 PetscCall(PetscMalloc1(n,&cnew)); 2968 for (i=0; i<m; i++) { 2969 PetscCall(MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2970 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2971 PetscCall(MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES)); 2972 PetscCall(MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2973 } 2974 PetscCall(PetscFree(cnew)); 2975 2976 (*B)->assembled = PETSC_FALSE; 2977 2978 #if defined(PETSC_HAVE_DEVICE) 2979 PetscCall(MatBindToCPU(*B,A->boundtocpu)); 2980 #endif 2981 PetscCall(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 2982 PetscCall(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY)); 2983 PetscCall(ISRestoreIndices(irowp,&row)); 2984 PetscCall(ISRestoreIndices(icolp,&col)); 2985 PetscCall(ISDestroy(&irowp)); 2986 PetscCall(ISDestroy(&icolp)); 2987 if (rowp == colp) { 2988 PetscCall(MatPropagateSymmetryOptions(A,*B)); 2989 } 2990 PetscFunctionReturn(0); 2991 } 2992 2993 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2994 { 2995 PetscFunctionBegin; 2996 /* If the two matrices have the same copy implementation, use fast copy. */ 2997 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2998 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2999 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 3000 const PetscScalar *aa; 3001 3002 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 3003 PetscCheck(a->i[A->rmap->n] == b->i[B->rmap->n],PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different %" PetscInt_FMT " != %" PetscInt_FMT,a->i[A->rmap->n],b->i[B->rmap->n]); 3004 PetscCall(PetscArraycpy(b->a,aa,a->i[A->rmap->n])); 3005 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 3006 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 3007 } else { 3008 PetscCall(MatCopy_Basic(A,B,str)); 3009 } 3010 PetscFunctionReturn(0); 3011 } 3012 3013 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 3014 { 3015 PetscFunctionBegin; 3016 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,NULL)); 3017 PetscFunctionReturn(0); 3018 } 3019 3020 PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 3021 { 3022 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3023 3024 PetscFunctionBegin; 3025 *array = a->a; 3026 PetscFunctionReturn(0); 3027 } 3028 3029 PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 3030 { 3031 PetscFunctionBegin; 3032 *array = NULL; 3033 PetscFunctionReturn(0); 3034 } 3035 3036 /* 3037 Computes the number of nonzeros per row needed for preallocation when X and Y 3038 have different nonzero structure. 3039 */ 3040 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 3041 { 3042 PetscInt i,j,k,nzx,nzy; 3043 3044 PetscFunctionBegin; 3045 /* Set the number of nonzeros in the new matrix */ 3046 for (i=0; i<m; i++) { 3047 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 3048 nzx = xi[i+1] - xi[i]; 3049 nzy = yi[i+1] - yi[i]; 3050 nnz[i] = 0; 3051 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 3052 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 3053 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 3054 nnz[i]++; 3055 } 3056 for (; k<nzy; k++) nnz[i]++; 3057 } 3058 PetscFunctionReturn(0); 3059 } 3060 3061 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 3062 { 3063 PetscInt m = Y->rmap->N; 3064 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 3065 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 3066 3067 PetscFunctionBegin; 3068 /* Set the number of nonzeros in the new matrix */ 3069 PetscCall(MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz)); 3070 PetscFunctionReturn(0); 3071 } 3072 3073 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 3074 { 3075 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 3076 3077 PetscFunctionBegin; 3078 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 3079 PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE; 3080 if (e) { 3081 PetscCall(PetscArraycmp(x->i,y->i,Y->rmap->n+1,&e)); 3082 if (e) { 3083 PetscCall(PetscArraycmp(x->j,y->j,y->nz,&e)); 3084 if (e) str = SAME_NONZERO_PATTERN; 3085 } 3086 } 3087 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatStructure is not SAME_NONZERO_PATTERN"); 3088 } 3089 if (str == SAME_NONZERO_PATTERN) { 3090 const PetscScalar *xa; 3091 PetscScalar *ya,alpha = a; 3092 PetscBLASInt one = 1,bnz; 3093 3094 PetscCall(PetscBLASIntCast(x->nz,&bnz)); 3095 PetscCall(MatSeqAIJGetArray(Y,&ya)); 3096 PetscCall(MatSeqAIJGetArrayRead(X,&xa)); 3097 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,xa,&one,ya,&one)); 3098 PetscCall(MatSeqAIJRestoreArrayRead(X,&xa)); 3099 PetscCall(MatSeqAIJRestoreArray(Y,&ya)); 3100 PetscCall(PetscLogFlops(2.0*bnz)); 3101 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 3102 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3103 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 3104 PetscCall(MatAXPY_Basic(Y,a,X,str)); 3105 } else { 3106 Mat B; 3107 PetscInt *nnz; 3108 PetscCall(PetscMalloc1(Y->rmap->N,&nnz)); 3109 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y),&B)); 3110 PetscCall(PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name)); 3111 PetscCall(MatSetLayouts(B,Y->rmap,Y->cmap)); 3112 PetscCall(MatSetType(B,((PetscObject)Y)->type_name)); 3113 PetscCall(MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz)); 3114 PetscCall(MatSeqAIJSetPreallocation(B,0,nnz)); 3115 PetscCall(MatAXPY_BasicWithPreallocation(B,Y,a,X,str)); 3116 PetscCall(MatHeaderMerge(Y,&B)); 3117 PetscCall(PetscFree(nnz)); 3118 } 3119 PetscFunctionReturn(0); 3120 } 3121 3122 PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 3123 { 3124 #if defined(PETSC_USE_COMPLEX) 3125 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3126 PetscInt i,nz; 3127 PetscScalar *a; 3128 3129 PetscFunctionBegin; 3130 nz = aij->nz; 3131 PetscCall(MatSeqAIJGetArray(mat,&a)); 3132 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 3133 PetscCall(MatSeqAIJRestoreArray(mat,&a)); 3134 #else 3135 PetscFunctionBegin; 3136 #endif 3137 PetscFunctionReturn(0); 3138 } 3139 3140 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3141 { 3142 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3143 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3144 PetscReal atmp; 3145 PetscScalar *x; 3146 const MatScalar *aa,*av; 3147 3148 PetscFunctionBegin; 3149 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3150 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3151 aa = av; 3152 ai = a->i; 3153 aj = a->j; 3154 3155 PetscCall(VecSet(v,0.0)); 3156 PetscCall(VecGetArrayWrite(v,&x)); 3157 PetscCall(VecGetLocalSize(v,&n)); 3158 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3159 for (i=0; i<m; i++) { 3160 ncols = ai[1] - ai[0]; ai++; 3161 for (j=0; j<ncols; j++) { 3162 atmp = PetscAbsScalar(*aa); 3163 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3164 aa++; aj++; 3165 } 3166 } 3167 PetscCall(VecRestoreArrayWrite(v,&x)); 3168 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3169 PetscFunctionReturn(0); 3170 } 3171 3172 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3173 { 3174 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3175 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3176 PetscScalar *x; 3177 const MatScalar *aa,*av; 3178 3179 PetscFunctionBegin; 3180 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3181 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3182 aa = av; 3183 ai = a->i; 3184 aj = a->j; 3185 3186 PetscCall(VecSet(v,0.0)); 3187 PetscCall(VecGetArrayWrite(v,&x)); 3188 PetscCall(VecGetLocalSize(v,&n)); 3189 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3190 for (i=0; i<m; i++) { 3191 ncols = ai[1] - ai[0]; ai++; 3192 if (ncols == A->cmap->n) { /* row is dense */ 3193 x[i] = *aa; if (idx) idx[i] = 0; 3194 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3195 x[i] = 0.0; 3196 if (idx) { 3197 for (j=0; j<ncols; j++) { /* find first implicit 0.0 in the row */ 3198 if (aj[j] > j) { 3199 idx[i] = j; 3200 break; 3201 } 3202 } 3203 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3204 if (j==ncols && j < A->cmap->n) idx[i] = j; 3205 } 3206 } 3207 for (j=0; j<ncols; j++) { 3208 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3209 aa++; aj++; 3210 } 3211 } 3212 PetscCall(VecRestoreArrayWrite(v,&x)); 3213 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3214 PetscFunctionReturn(0); 3215 } 3216 3217 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3218 { 3219 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3220 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3221 PetscScalar *x; 3222 const MatScalar *aa,*av; 3223 3224 PetscFunctionBegin; 3225 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3226 aa = av; 3227 ai = a->i; 3228 aj = a->j; 3229 3230 PetscCall(VecSet(v,0.0)); 3231 PetscCall(VecGetArrayWrite(v,&x)); 3232 PetscCall(VecGetLocalSize(v,&n)); 3233 PetscCheck(n == m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %" PetscInt_FMT " vs. %" PetscInt_FMT " rows", m, n); 3234 for (i=0; i<m; i++) { 3235 ncols = ai[1] - ai[0]; ai++; 3236 if (ncols == A->cmap->n) { /* row is dense */ 3237 x[i] = *aa; if (idx) idx[i] = 0; 3238 } else { /* row is sparse so already KNOW minimum is 0.0 or higher */ 3239 x[i] = 0.0; 3240 if (idx) { /* find first implicit 0.0 in the row */ 3241 for (j=0; j<ncols; j++) { 3242 if (aj[j] > j) { 3243 idx[i] = j; 3244 break; 3245 } 3246 } 3247 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3248 if (j==ncols && j < A->cmap->n) idx[i] = j; 3249 } 3250 } 3251 for (j=0; j<ncols; j++) { 3252 if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3253 aa++; aj++; 3254 } 3255 } 3256 PetscCall(VecRestoreArrayWrite(v,&x)); 3257 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3258 PetscFunctionReturn(0); 3259 } 3260 3261 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3262 { 3263 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3264 PetscInt i,j,m = A->rmap->n,ncols,n; 3265 const PetscInt *ai,*aj; 3266 PetscScalar *x; 3267 const MatScalar *aa,*av; 3268 3269 PetscFunctionBegin; 3270 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3271 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3272 aa = av; 3273 ai = a->i; 3274 aj = a->j; 3275 3276 PetscCall(VecSet(v,0.0)); 3277 PetscCall(VecGetArrayWrite(v,&x)); 3278 PetscCall(VecGetLocalSize(v,&n)); 3279 PetscCheck(n == m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3280 for (i=0; i<m; i++) { 3281 ncols = ai[1] - ai[0]; ai++; 3282 if (ncols == A->cmap->n) { /* row is dense */ 3283 x[i] = *aa; if (idx) idx[i] = 0; 3284 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3285 x[i] = 0.0; 3286 if (idx) { /* find first implicit 0.0 in the row */ 3287 for (j=0; j<ncols; j++) { 3288 if (aj[j] > j) { 3289 idx[i] = j; 3290 break; 3291 } 3292 } 3293 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3294 if (j==ncols && j < A->cmap->n) idx[i] = j; 3295 } 3296 } 3297 for (j=0; j<ncols; j++) { 3298 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3299 aa++; aj++; 3300 } 3301 } 3302 PetscCall(VecRestoreArrayWrite(v,&x)); 3303 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3304 PetscFunctionReturn(0); 3305 } 3306 3307 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3308 { 3309 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3310 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3311 MatScalar *diag,work[25],*v_work; 3312 const PetscReal shift = 0.0; 3313 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3314 3315 PetscFunctionBegin; 3316 allowzeropivot = PetscNot(A->erroriffailure); 3317 if (a->ibdiagvalid) { 3318 if (values) *values = a->ibdiag; 3319 PetscFunctionReturn(0); 3320 } 3321 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 3322 if (!a->ibdiag) { 3323 PetscCall(PetscMalloc1(bs2*mbs,&a->ibdiag)); 3324 PetscCall(PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar))); 3325 } 3326 diag = a->ibdiag; 3327 if (values) *values = a->ibdiag; 3328 /* factor and invert each block */ 3329 switch (bs) { 3330 case 1: 3331 for (i=0; i<mbs; i++) { 3332 PetscCall(MatGetValues(A,1,&i,1,&i,diag+i)); 3333 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3334 if (allowzeropivot) { 3335 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3336 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3337 A->factorerror_zeropivot_row = i; 3338 PetscCall(PetscInfo(A,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON)); 3339 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON); 3340 } 3341 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3342 } 3343 break; 3344 case 2: 3345 for (i=0; i<mbs; i++) { 3346 ij[0] = 2*i; ij[1] = 2*i + 1; 3347 PetscCall(MatGetValues(A,2,ij,2,ij,diag)); 3348 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 3349 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3350 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 3351 diag += 4; 3352 } 3353 break; 3354 case 3: 3355 for (i=0; i<mbs; i++) { 3356 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3357 PetscCall(MatGetValues(A,3,ij,3,ij,diag)); 3358 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 3359 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3360 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 3361 diag += 9; 3362 } 3363 break; 3364 case 4: 3365 for (i=0; i<mbs; i++) { 3366 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3367 PetscCall(MatGetValues(A,4,ij,4,ij,diag)); 3368 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 3369 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3370 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 3371 diag += 16; 3372 } 3373 break; 3374 case 5: 3375 for (i=0; i<mbs; i++) { 3376 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3377 PetscCall(MatGetValues(A,5,ij,5,ij,diag)); 3378 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 3379 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3380 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 3381 diag += 25; 3382 } 3383 break; 3384 case 6: 3385 for (i=0; i<mbs; i++) { 3386 ij[0] = 6*i; ij[1] = 6*i + 1; ij[2] = 6*i + 2; ij[3] = 6*i + 3; ij[4] = 6*i + 4; ij[5] = 6*i + 5; 3387 PetscCall(MatGetValues(A,6,ij,6,ij,diag)); 3388 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 3389 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3390 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 3391 diag += 36; 3392 } 3393 break; 3394 case 7: 3395 for (i=0; i<mbs; i++) { 3396 ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6; 3397 PetscCall(MatGetValues(A,7,ij,7,ij,diag)); 3398 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 3399 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3400 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 3401 diag += 49; 3402 } 3403 break; 3404 default: 3405 PetscCall(PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ)); 3406 for (i=0; i<mbs; i++) { 3407 for (j=0; j<bs; j++) { 3408 IJ[j] = bs*i + j; 3409 } 3410 PetscCall(MatGetValues(A,bs,IJ,bs,IJ,diag)); 3411 PetscCall(PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 3412 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3413 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bs)); 3414 diag += bs2; 3415 } 3416 PetscCall(PetscFree3(v_work,v_pivots,IJ)); 3417 } 3418 a->ibdiagvalid = PETSC_TRUE; 3419 PetscFunctionReturn(0); 3420 } 3421 3422 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3423 { 3424 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3425 PetscScalar a,*aa; 3426 PetscInt m,n,i,j,col; 3427 3428 PetscFunctionBegin; 3429 if (!x->assembled) { 3430 PetscCall(MatGetSize(x,&m,&n)); 3431 for (i=0; i<m; i++) { 3432 for (j=0; j<aij->imax[i]; j++) { 3433 PetscCall(PetscRandomGetValue(rctx,&a)); 3434 col = (PetscInt)(n*PetscRealPart(a)); 3435 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3436 } 3437 } 3438 } else { 3439 PetscCall(MatSeqAIJGetArrayWrite(x,&aa)); 3440 for (i=0; i<aij->nz; i++) PetscCall(PetscRandomGetValue(rctx,aa+i)); 3441 PetscCall(MatSeqAIJRestoreArrayWrite(x,&aa)); 3442 } 3443 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3444 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3445 PetscFunctionReturn(0); 3446 } 3447 3448 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3449 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx) 3450 { 3451 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3452 PetscScalar a; 3453 PetscInt m,n,i,j,col,nskip; 3454 3455 PetscFunctionBegin; 3456 nskip = high - low; 3457 PetscCall(MatGetSize(x,&m,&n)); 3458 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3459 for (i=0; i<m; i++) { 3460 for (j=0; j<aij->imax[i]; j++) { 3461 PetscCall(PetscRandomGetValue(rctx,&a)); 3462 col = (PetscInt)(n*PetscRealPart(a)); 3463 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3464 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3465 } 3466 } 3467 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3468 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3469 PetscFunctionReturn(0); 3470 } 3471 3472 /* -------------------------------------------------------------------*/ 3473 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3474 MatGetRow_SeqAIJ, 3475 MatRestoreRow_SeqAIJ, 3476 MatMult_SeqAIJ, 3477 /* 4*/ MatMultAdd_SeqAIJ, 3478 MatMultTranspose_SeqAIJ, 3479 MatMultTransposeAdd_SeqAIJ, 3480 NULL, 3481 NULL, 3482 NULL, 3483 /* 10*/ NULL, 3484 MatLUFactor_SeqAIJ, 3485 NULL, 3486 MatSOR_SeqAIJ, 3487 MatTranspose_SeqAIJ, 3488 /*1 5*/ MatGetInfo_SeqAIJ, 3489 MatEqual_SeqAIJ, 3490 MatGetDiagonal_SeqAIJ, 3491 MatDiagonalScale_SeqAIJ, 3492 MatNorm_SeqAIJ, 3493 /* 20*/ NULL, 3494 MatAssemblyEnd_SeqAIJ, 3495 MatSetOption_SeqAIJ, 3496 MatZeroEntries_SeqAIJ, 3497 /* 24*/ MatZeroRows_SeqAIJ, 3498 NULL, 3499 NULL, 3500 NULL, 3501 NULL, 3502 /* 29*/ MatSetUp_SeqAIJ, 3503 NULL, 3504 NULL, 3505 NULL, 3506 NULL, 3507 /* 34*/ MatDuplicate_SeqAIJ, 3508 NULL, 3509 NULL, 3510 MatILUFactor_SeqAIJ, 3511 NULL, 3512 /* 39*/ MatAXPY_SeqAIJ, 3513 MatCreateSubMatrices_SeqAIJ, 3514 MatIncreaseOverlap_SeqAIJ, 3515 MatGetValues_SeqAIJ, 3516 MatCopy_SeqAIJ, 3517 /* 44*/ MatGetRowMax_SeqAIJ, 3518 MatScale_SeqAIJ, 3519 MatShift_SeqAIJ, 3520 MatDiagonalSet_SeqAIJ, 3521 MatZeroRowsColumns_SeqAIJ, 3522 /* 49*/ MatSetRandom_SeqAIJ, 3523 MatGetRowIJ_SeqAIJ, 3524 MatRestoreRowIJ_SeqAIJ, 3525 MatGetColumnIJ_SeqAIJ, 3526 MatRestoreColumnIJ_SeqAIJ, 3527 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3528 NULL, 3529 NULL, 3530 MatPermute_SeqAIJ, 3531 NULL, 3532 /* 59*/ NULL, 3533 MatDestroy_SeqAIJ, 3534 MatView_SeqAIJ, 3535 NULL, 3536 NULL, 3537 /* 64*/ NULL, 3538 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3539 NULL, 3540 NULL, 3541 NULL, 3542 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3543 MatGetRowMinAbs_SeqAIJ, 3544 NULL, 3545 NULL, 3546 NULL, 3547 /* 74*/ NULL, 3548 MatFDColoringApply_AIJ, 3549 NULL, 3550 NULL, 3551 NULL, 3552 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3553 NULL, 3554 NULL, 3555 NULL, 3556 MatLoad_SeqAIJ, 3557 /* 84*/ MatIsSymmetric_SeqAIJ, 3558 MatIsHermitian_SeqAIJ, 3559 NULL, 3560 NULL, 3561 NULL, 3562 /* 89*/ NULL, 3563 NULL, 3564 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3565 NULL, 3566 NULL, 3567 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3568 NULL, 3569 NULL, 3570 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3571 NULL, 3572 /* 99*/ MatProductSetFromOptions_SeqAIJ, 3573 NULL, 3574 NULL, 3575 MatConjugate_SeqAIJ, 3576 NULL, 3577 /*104*/ MatSetValuesRow_SeqAIJ, 3578 MatRealPart_SeqAIJ, 3579 MatImaginaryPart_SeqAIJ, 3580 NULL, 3581 NULL, 3582 /*109*/ MatMatSolve_SeqAIJ, 3583 NULL, 3584 MatGetRowMin_SeqAIJ, 3585 NULL, 3586 MatMissingDiagonal_SeqAIJ, 3587 /*114*/ NULL, 3588 NULL, 3589 NULL, 3590 NULL, 3591 NULL, 3592 /*119*/ NULL, 3593 NULL, 3594 NULL, 3595 NULL, 3596 MatGetMultiProcBlock_SeqAIJ, 3597 /*124*/ MatFindNonzeroRows_SeqAIJ, 3598 MatGetColumnReductions_SeqAIJ, 3599 MatInvertBlockDiagonal_SeqAIJ, 3600 MatInvertVariableBlockDiagonal_SeqAIJ, 3601 NULL, 3602 /*129*/ NULL, 3603 NULL, 3604 NULL, 3605 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3606 MatTransposeColoringCreate_SeqAIJ, 3607 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3608 MatTransColoringApplyDenToSp_SeqAIJ, 3609 NULL, 3610 NULL, 3611 MatRARtNumeric_SeqAIJ_SeqAIJ, 3612 /*139*/NULL, 3613 NULL, 3614 NULL, 3615 MatFDColoringSetUp_SeqXAIJ, 3616 MatFindOffBlockDiagonalEntries_SeqAIJ, 3617 MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3618 /*145*/MatDestroySubMatrices_SeqAIJ, 3619 NULL, 3620 NULL 3621 }; 3622 3623 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3624 { 3625 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3626 PetscInt i,nz,n; 3627 3628 PetscFunctionBegin; 3629 nz = aij->maxnz; 3630 n = mat->rmap->n; 3631 for (i=0; i<nz; i++) { 3632 aij->j[i] = indices[i]; 3633 } 3634 aij->nz = nz; 3635 for (i=0; i<n; i++) { 3636 aij->ilen[i] = aij->imax[i]; 3637 } 3638 PetscFunctionReturn(0); 3639 } 3640 3641 /* 3642 * Given a sparse matrix with global column indices, compact it by using a local column space. 3643 * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3644 */ 3645 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3646 { 3647 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3648 PetscTable gid1_lid1; 3649 PetscTablePosition tpos; 3650 PetscInt gid,lid,i,ec,nz = aij->nz; 3651 PetscInt *garray,*jj = aij->j; 3652 3653 PetscFunctionBegin; 3654 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3655 PetscValidPointer(mapping,2); 3656 /* use a table */ 3657 PetscCall(PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1)); 3658 ec = 0; 3659 for (i=0; i<nz; i++) { 3660 PetscInt data,gid1 = jj[i] + 1; 3661 PetscCall(PetscTableFind(gid1_lid1,gid1,&data)); 3662 if (!data) { 3663 /* one based table */ 3664 PetscCall(PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES)); 3665 } 3666 } 3667 /* form array of columns we need */ 3668 PetscCall(PetscMalloc1(ec,&garray)); 3669 PetscCall(PetscTableGetHeadPosition(gid1_lid1,&tpos)); 3670 while (tpos) { 3671 PetscCall(PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid)); 3672 gid--; 3673 lid--; 3674 garray[lid] = gid; 3675 } 3676 PetscCall(PetscSortInt(ec,garray)); /* sort, and rebuild */ 3677 PetscCall(PetscTableRemoveAll(gid1_lid1)); 3678 for (i=0; i<ec; i++) { 3679 PetscCall(PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES)); 3680 } 3681 /* compact out the extra columns in B */ 3682 for (i=0; i<nz; i++) { 3683 PetscInt gid1 = jj[i] + 1; 3684 PetscCall(PetscTableFind(gid1_lid1,gid1,&lid)); 3685 lid--; 3686 jj[i] = lid; 3687 } 3688 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3689 PetscCall(PetscTableDestroy(&gid1_lid1)); 3690 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat),ec,ec,1,&mat->cmap)); 3691 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping)); 3692 PetscCall(ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH)); 3693 PetscFunctionReturn(0); 3694 } 3695 3696 /*@ 3697 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3698 in the matrix. 3699 3700 Input Parameters: 3701 + mat - the SeqAIJ matrix 3702 - indices - the column indices 3703 3704 Level: advanced 3705 3706 Notes: 3707 This can be called if you have precomputed the nonzero structure of the 3708 matrix and want to provide it to the matrix object to improve the performance 3709 of the MatSetValues() operation. 3710 3711 You MUST have set the correct numbers of nonzeros per row in the call to 3712 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3713 3714 MUST be called before any calls to MatSetValues(); 3715 3716 The indices should start with zero, not one. 3717 3718 @*/ 3719 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3720 { 3721 PetscFunctionBegin; 3722 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3723 PetscValidIntPointer(indices,2); 3724 PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices)); 3725 PetscFunctionReturn(0); 3726 } 3727 3728 /* ----------------------------------------------------------------------------------------*/ 3729 3730 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3731 { 3732 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3733 size_t nz = aij->i[mat->rmap->n]; 3734 3735 PetscFunctionBegin; 3736 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3737 3738 /* allocate space for values if not already there */ 3739 if (!aij->saved_values) { 3740 PetscCall(PetscMalloc1(nz+1,&aij->saved_values)); 3741 PetscCall(PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar))); 3742 } 3743 3744 /* copy values over */ 3745 PetscCall(PetscArraycpy(aij->saved_values,aij->a,nz)); 3746 PetscFunctionReturn(0); 3747 } 3748 3749 /*@ 3750 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3751 example, reuse of the linear part of a Jacobian, while recomputing the 3752 nonlinear portion. 3753 3754 Collect on Mat 3755 3756 Input Parameters: 3757 . mat - the matrix (currently only AIJ matrices support this option) 3758 3759 Level: advanced 3760 3761 Common Usage, with SNESSolve(): 3762 $ Create Jacobian matrix 3763 $ Set linear terms into matrix 3764 $ Apply boundary conditions to matrix, at this time matrix must have 3765 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3766 $ boundary conditions again will not change the nonzero structure 3767 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3768 $ ierr = MatStoreValues(mat); 3769 $ Call SNESSetJacobian() with matrix 3770 $ In your Jacobian routine 3771 $ ierr = MatRetrieveValues(mat); 3772 $ Set nonlinear terms in matrix 3773 3774 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3775 $ // build linear portion of Jacobian 3776 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3777 $ ierr = MatStoreValues(mat); 3778 $ loop over nonlinear iterations 3779 $ ierr = MatRetrieveValues(mat); 3780 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3781 $ // call MatAssemblyBegin/End() on matrix 3782 $ Solve linear system with Jacobian 3783 $ endloop 3784 3785 Notes: 3786 Matrix must already be assemblied before calling this routine 3787 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3788 calling this routine. 3789 3790 When this is called multiple times it overwrites the previous set of stored values 3791 and does not allocated additional space. 3792 3793 .seealso: MatRetrieveValues() 3794 3795 @*/ 3796 PetscErrorCode MatStoreValues(Mat mat) 3797 { 3798 PetscFunctionBegin; 3799 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3800 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3801 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3802 PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat)); 3803 PetscFunctionReturn(0); 3804 } 3805 3806 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3807 { 3808 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3809 PetscInt nz = aij->i[mat->rmap->n]; 3810 3811 PetscFunctionBegin; 3812 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3813 PetscCheck(aij->saved_values,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3814 /* copy values over */ 3815 PetscCall(PetscArraycpy(aij->a,aij->saved_values,nz)); 3816 PetscFunctionReturn(0); 3817 } 3818 3819 /*@ 3820 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3821 example, reuse of the linear part of a Jacobian, while recomputing the 3822 nonlinear portion. 3823 3824 Collect on Mat 3825 3826 Input Parameters: 3827 . mat - the matrix (currently only AIJ matrices support this option) 3828 3829 Level: advanced 3830 3831 .seealso: MatStoreValues() 3832 3833 @*/ 3834 PetscErrorCode MatRetrieveValues(Mat mat) 3835 { 3836 PetscFunctionBegin; 3837 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3838 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3839 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3840 PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat)); 3841 PetscFunctionReturn(0); 3842 } 3843 3844 /* --------------------------------------------------------------------------------*/ 3845 /*@C 3846 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3847 (the default parallel PETSc format). For good matrix assembly performance 3848 the user should preallocate the matrix storage by setting the parameter nz 3849 (or the array nnz). By setting these parameters accurately, performance 3850 during matrix assembly can be increased by more than a factor of 50. 3851 3852 Collective 3853 3854 Input Parameters: 3855 + comm - MPI communicator, set to PETSC_COMM_SELF 3856 . m - number of rows 3857 . n - number of columns 3858 . nz - number of nonzeros per row (same for all rows) 3859 - nnz - array containing the number of nonzeros in the various rows 3860 (possibly different for each row) or NULL 3861 3862 Output Parameter: 3863 . A - the matrix 3864 3865 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3866 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3867 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3868 3869 Notes: 3870 If nnz is given then nz is ignored 3871 3872 The AIJ format (also called the Yale sparse matrix format or 3873 compressed row storage), is fully compatible with standard Fortran 77 3874 storage. That is, the stored row and column indices can begin at 3875 either one (as in Fortran) or zero. See the users' manual for details. 3876 3877 Specify the preallocated storage with either nz or nnz (not both). 3878 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3879 allocation. For large problems you MUST preallocate memory or you 3880 will get TERRIBLE performance, see the users' manual chapter on matrices. 3881 3882 By default, this format uses inodes (identical nodes) when possible, to 3883 improve numerical efficiency of matrix-vector products and solves. We 3884 search for consecutive rows with the same nonzero structure, thereby 3885 reusing matrix information to achieve increased efficiency. 3886 3887 Options Database Keys: 3888 + -mat_no_inode - Do not use inodes 3889 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3890 3891 Level: intermediate 3892 3893 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3894 3895 @*/ 3896 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3897 { 3898 PetscFunctionBegin; 3899 PetscCall(MatCreate(comm,A)); 3900 PetscCall(MatSetSizes(*A,m,n,m,n)); 3901 PetscCall(MatSetType(*A,MATSEQAIJ)); 3902 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz)); 3903 PetscFunctionReturn(0); 3904 } 3905 3906 /*@C 3907 MatSeqAIJSetPreallocation - For good matrix assembly performance 3908 the user should preallocate the matrix storage by setting the parameter nz 3909 (or the array nnz). By setting these parameters accurately, performance 3910 during matrix assembly can be increased by more than a factor of 50. 3911 3912 Collective 3913 3914 Input Parameters: 3915 + B - The matrix 3916 . nz - number of nonzeros per row (same for all rows) 3917 - nnz - array containing the number of nonzeros in the various rows 3918 (possibly different for each row) or NULL 3919 3920 Notes: 3921 If nnz is given then nz is ignored 3922 3923 The AIJ format (also called the Yale sparse matrix format or 3924 compressed row storage), is fully compatible with standard Fortran 77 3925 storage. That is, the stored row and column indices can begin at 3926 either one (as in Fortran) or zero. See the users' manual for details. 3927 3928 Specify the preallocated storage with either nz or nnz (not both). 3929 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3930 allocation. For large problems you MUST preallocate memory or you 3931 will get TERRIBLE performance, see the users' manual chapter on matrices. 3932 3933 You can call MatGetInfo() to get information on how effective the preallocation was; 3934 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3935 You can also run with the option -info and look for messages with the string 3936 malloc in them to see if additional memory allocation was needed. 3937 3938 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3939 entries or columns indices 3940 3941 By default, this format uses inodes (identical nodes) when possible, to 3942 improve numerical efficiency of matrix-vector products and solves. We 3943 search for consecutive rows with the same nonzero structure, thereby 3944 reusing matrix information to achieve increased efficiency. 3945 3946 Options Database Keys: 3947 + -mat_no_inode - Do not use inodes 3948 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3949 3950 Level: intermediate 3951 3952 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo(), 3953 MatSeqAIJSetTotalPreallocation() 3954 3955 @*/ 3956 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3957 { 3958 PetscFunctionBegin; 3959 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3960 PetscValidType(B,1); 3961 PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz)); 3962 PetscFunctionReturn(0); 3963 } 3964 3965 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3966 { 3967 Mat_SeqAIJ *b; 3968 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3969 PetscInt i; 3970 3971 PetscFunctionBegin; 3972 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3973 if (nz == MAT_SKIP_ALLOCATION) { 3974 skipallocation = PETSC_TRUE; 3975 nz = 0; 3976 } 3977 PetscCall(PetscLayoutSetUp(B->rmap)); 3978 PetscCall(PetscLayoutSetUp(B->cmap)); 3979 3980 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3981 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %" PetscInt_FMT,nz); 3982 if (PetscUnlikelyDebug(nnz)) { 3983 for (i=0; i<B->rmap->n; i++) { 3984 PetscCheck(nnz[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT,i,nnz[i]); 3985 PetscCheck(nnz[i] <= B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " rowlength %" PetscInt_FMT,i,nnz[i],B->cmap->n); 3986 } 3987 } 3988 3989 B->preallocated = PETSC_TRUE; 3990 3991 b = (Mat_SeqAIJ*)B->data; 3992 3993 if (!skipallocation) { 3994 if (!b->imax) { 3995 PetscCall(PetscMalloc1(B->rmap->n,&b->imax)); 3996 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 3997 } 3998 if (!b->ilen) { 3999 /* b->ilen will count nonzeros in each row so far. */ 4000 PetscCall(PetscCalloc1(B->rmap->n,&b->ilen)); 4001 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4002 } else { 4003 PetscCall(PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt))); 4004 } 4005 if (!b->ipre) { 4006 PetscCall(PetscMalloc1(B->rmap->n,&b->ipre)); 4007 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4008 } 4009 if (!nnz) { 4010 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 4011 else if (nz < 0) nz = 1; 4012 nz = PetscMin(nz,B->cmap->n); 4013 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 4014 nz = nz*B->rmap->n; 4015 } else { 4016 PetscInt64 nz64 = 0; 4017 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 4018 PetscCall(PetscIntCast(nz64,&nz)); 4019 } 4020 4021 /* allocate the matrix space */ 4022 /* FIXME: should B's old memory be unlogged? */ 4023 PetscCall(MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i)); 4024 if (B->structure_only) { 4025 PetscCall(PetscMalloc1(nz,&b->j)); 4026 PetscCall(PetscMalloc1(B->rmap->n+1,&b->i)); 4027 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt))); 4028 } else { 4029 PetscCall(PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i)); 4030 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)))); 4031 } 4032 b->i[0] = 0; 4033 for (i=1; i<B->rmap->n+1; i++) { 4034 b->i[i] = b->i[i-1] + b->imax[i-1]; 4035 } 4036 if (B->structure_only) { 4037 b->singlemalloc = PETSC_FALSE; 4038 b->free_a = PETSC_FALSE; 4039 } else { 4040 b->singlemalloc = PETSC_TRUE; 4041 b->free_a = PETSC_TRUE; 4042 } 4043 b->free_ij = PETSC_TRUE; 4044 } else { 4045 b->free_a = PETSC_FALSE; 4046 b->free_ij = PETSC_FALSE; 4047 } 4048 4049 if (b->ipre && nnz != b->ipre && b->imax) { 4050 /* reserve user-requested sparsity */ 4051 PetscCall(PetscArraycpy(b->ipre,b->imax,B->rmap->n)); 4052 } 4053 4054 b->nz = 0; 4055 b->maxnz = nz; 4056 B->info.nz_unneeded = (double)b->maxnz; 4057 if (realalloc) { 4058 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 4059 } 4060 B->was_assembled = PETSC_FALSE; 4061 B->assembled = PETSC_FALSE; 4062 PetscFunctionReturn(0); 4063 } 4064 4065 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4066 { 4067 Mat_SeqAIJ *a; 4068 PetscInt i; 4069 4070 PetscFunctionBegin; 4071 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4072 4073 /* Check local size. If zero, then return */ 4074 if (!A->rmap->n) PetscFunctionReturn(0); 4075 4076 a = (Mat_SeqAIJ*)A->data; 4077 /* if no saved info, we error out */ 4078 PetscCheck(a->ipre,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"No saved preallocation info "); 4079 4080 PetscCheckFalse(!a->i || !a->j || !a->a || !a->imax || !a->ilen,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Memory info is incomplete, and can not reset preallocation "); 4081 4082 PetscCall(PetscArraycpy(a->imax,a->ipre,A->rmap->n)); 4083 PetscCall(PetscArrayzero(a->ilen,A->rmap->n)); 4084 a->i[0] = 0; 4085 for (i=1; i<A->rmap->n+1; i++) { 4086 a->i[i] = a->i[i-1] + a->imax[i-1]; 4087 } 4088 A->preallocated = PETSC_TRUE; 4089 a->nz = 0; 4090 a->maxnz = a->i[A->rmap->n]; 4091 A->info.nz_unneeded = (double)a->maxnz; 4092 A->was_assembled = PETSC_FALSE; 4093 A->assembled = PETSC_FALSE; 4094 PetscFunctionReturn(0); 4095 } 4096 4097 /*@ 4098 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 4099 4100 Input Parameters: 4101 + B - the matrix 4102 . i - the indices into j for the start of each row (starts with zero) 4103 . j - the column indices for each row (starts with zero) these must be sorted for each row 4104 - v - optional values in the matrix 4105 4106 Level: developer 4107 4108 Notes: 4109 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 4110 4111 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4112 structure will be the union of all the previous nonzero structures. 4113 4114 Developer Notes: 4115 An optimization could be added to the implementation where it checks if the i, and j are identical to the current i and j and 4116 then just copies the v values directly with PetscMemcpy(). 4117 4118 This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them. 4119 4120 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ, MatResetPreallocation() 4121 @*/ 4122 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 4123 { 4124 PetscFunctionBegin; 4125 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 4126 PetscValidType(B,1); 4127 PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v)); 4128 PetscFunctionReturn(0); 4129 } 4130 4131 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 4132 { 4133 PetscInt i; 4134 PetscInt m,n; 4135 PetscInt nz; 4136 PetscInt *nnz; 4137 4138 PetscFunctionBegin; 4139 PetscCheckFalse(Ii[0],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4140 4141 PetscCall(PetscLayoutSetUp(B->rmap)); 4142 PetscCall(PetscLayoutSetUp(B->cmap)); 4143 4144 PetscCall(MatGetSize(B, &m, &n)); 4145 PetscCall(PetscMalloc1(m+1, &nnz)); 4146 for (i = 0; i < m; i++) { 4147 nz = Ii[i+1]- Ii[i]; 4148 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4149 nnz[i] = nz; 4150 } 4151 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4152 PetscCall(PetscFree(nnz)); 4153 4154 for (i = 0; i < m; i++) { 4155 PetscCall(MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES)); 4156 } 4157 4158 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 4159 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 4160 4161 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 4162 PetscFunctionReturn(0); 4163 } 4164 4165 /*@ 4166 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4167 4168 Input Parameters: 4169 + A - left-hand side matrix 4170 . B - right-hand side matrix 4171 - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 4172 4173 Output Parameter: 4174 . C - Kronecker product of A and B 4175 4176 Level: intermediate 4177 4178 Notes: 4179 MAT_REUSE_MATRIX can only be used when the nonzero structure of the product matrix has not changed from that last call to MatSeqAIJKron(). 4180 4181 .seealso: MatCreateSeqAIJ(), MATSEQAIJ, MATKAIJ, MatReuse 4182 @*/ 4183 PetscErrorCode MatSeqAIJKron(Mat A,Mat B,MatReuse reuse,Mat *C) 4184 { 4185 PetscFunctionBegin; 4186 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4187 PetscValidType(A,1); 4188 PetscValidHeaderSpecific(B,MAT_CLASSID,2); 4189 PetscValidType(B,2); 4190 PetscValidPointer(C,4); 4191 if (reuse == MAT_REUSE_MATRIX) { 4192 PetscValidHeaderSpecific(*C,MAT_CLASSID,4); 4193 PetscValidType(*C,4); 4194 } 4195 PetscTryMethod(A,"MatSeqAIJKron_C",(Mat,Mat,MatReuse,Mat*),(A,B,reuse,C)); 4196 PetscFunctionReturn(0); 4197 } 4198 4199 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A,Mat B,MatReuse reuse,Mat *C) 4200 { 4201 Mat newmat; 4202 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4203 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 4204 PetscScalar *v; 4205 const PetscScalar *aa,*ba; 4206 PetscInt *i,*j,m,n,p,q,nnz = 0,am = A->rmap->n,bm = B->rmap->n,an = A->cmap->n, bn = B->cmap->n; 4207 PetscBool flg; 4208 4209 PetscFunctionBegin; 4210 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4211 PetscCheck(A->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4212 PetscCheck(!B->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4213 PetscCheck(B->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4214 PetscCall(PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&flg)); 4215 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatType %s",((PetscObject)B)->type_name); 4216 PetscCheckFalse(reuse != MAT_INITIAL_MATRIX && reuse != MAT_REUSE_MATRIX,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatReuse %d",(int)reuse); 4217 if (reuse == MAT_INITIAL_MATRIX) { 4218 PetscCall(PetscMalloc2(am*bm+1,&i,a->i[am]*b->i[bm],&j)); 4219 PetscCall(MatCreate(PETSC_COMM_SELF,&newmat)); 4220 PetscCall(MatSetSizes(newmat,am*bm,an*bn,am*bm,an*bn)); 4221 PetscCall(MatSetType(newmat,MATAIJ)); 4222 i[0] = 0; 4223 for (m = 0; m < am; ++m) { 4224 for (p = 0; p < bm; ++p) { 4225 i[m*bm + p + 1] = i[m*bm + p] + (a->i[m+1] - a->i[m]) * (b->i[p+1] - b->i[p]); 4226 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4227 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4228 j[nnz++] = a->j[n]*bn + b->j[q]; 4229 } 4230 } 4231 } 4232 } 4233 PetscCall(MatSeqAIJSetPreallocationCSR(newmat,i,j,NULL)); 4234 *C = newmat; 4235 PetscCall(PetscFree2(i,j)); 4236 nnz = 0; 4237 } 4238 PetscCall(MatSeqAIJGetArray(*C,&v)); 4239 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4240 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 4241 for (m = 0; m < am; ++m) { 4242 for (p = 0; p < bm; ++p) { 4243 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4244 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4245 v[nnz++] = aa[n] * ba[q]; 4246 } 4247 } 4248 } 4249 } 4250 PetscCall(MatSeqAIJRestoreArray(*C,&v)); 4251 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 4252 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 4253 PetscFunctionReturn(0); 4254 } 4255 4256 #include <../src/mat/impls/dense/seq/dense.h> 4257 #include <petsc/private/kernels/petscaxpy.h> 4258 4259 /* 4260 Computes (B'*A')' since computing B*A directly is untenable 4261 4262 n p p 4263 [ ] [ ] [ ] 4264 m [ A ] * n [ B ] = m [ C ] 4265 [ ] [ ] [ ] 4266 4267 */ 4268 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 4269 { 4270 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 4271 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 4272 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 4273 PetscInt i,j,n,m,q,p; 4274 const PetscInt *ii,*idx; 4275 const PetscScalar *b,*a,*a_q; 4276 PetscScalar *c,*c_q; 4277 PetscInt clda = sub_c->lda; 4278 PetscInt alda = sub_a->lda; 4279 4280 PetscFunctionBegin; 4281 m = A->rmap->n; 4282 n = A->cmap->n; 4283 p = B->cmap->n; 4284 a = sub_a->v; 4285 b = sub_b->a; 4286 c = sub_c->v; 4287 if (clda == m) { 4288 PetscCall(PetscArrayzero(c,m*p)); 4289 } else { 4290 for (j=0;j<p;j++) 4291 for (i=0;i<m;i++) 4292 c[j*clda + i] = 0.0; 4293 } 4294 ii = sub_b->i; 4295 idx = sub_b->j; 4296 for (i=0; i<n; i++) { 4297 q = ii[i+1] - ii[i]; 4298 while (q-->0) { 4299 c_q = c + clda*(*idx); 4300 a_q = a + alda*i; 4301 PetscKernelAXPY(c_q,*b,a_q,m); 4302 idx++; 4303 b++; 4304 } 4305 } 4306 PetscFunctionReturn(0); 4307 } 4308 4309 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C) 4310 { 4311 PetscInt m=A->rmap->n,n=B->cmap->n; 4312 PetscBool cisdense; 4313 4314 PetscFunctionBegin; 4315 PetscCheck(A->cmap->n == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %" PetscInt_FMT " != B->rmap->n %" PetscInt_FMT,A->cmap->n,B->rmap->n); 4316 PetscCall(MatSetSizes(C,m,n,m,n)); 4317 PetscCall(MatSetBlockSizesFromMats(C,A,B)); 4318 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 4319 if (!cisdense) { 4320 PetscCall(MatSetType(C,MATDENSE)); 4321 } 4322 PetscCall(MatSetUp(C)); 4323 4324 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4325 PetscFunctionReturn(0); 4326 } 4327 4328 /* ----------------------------------------------------------------*/ 4329 /*MC 4330 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4331 based on compressed sparse row format. 4332 4333 Options Database Keys: 4334 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4335 4336 Level: beginner 4337 4338 Notes: 4339 MatSetValues() may be called for this matrix type with a NULL argument for the numerical values, 4340 in this case the values associated with the rows and columns one passes in are set to zero 4341 in the matrix 4342 4343 MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no 4344 space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored 4345 4346 Developer Notes: 4347 It would be nice if all matrix formats supported passing NULL in for the numerical values 4348 4349 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType, MATSELL, MATSEQSELL, MATMPISELL 4350 M*/ 4351 4352 /*MC 4353 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4354 4355 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4356 and MATMPIAIJ otherwise. As a result, for single process communicators, 4357 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4358 for communicators controlling multiple processes. It is recommended that you call both of 4359 the above preallocation routines for simplicity. 4360 4361 Options Database Keys: 4362 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4363 4364 Developer Notes: 4365 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4366 enough exist. 4367 4368 Level: beginner 4369 4370 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ, MATMPIAIJ, MATSELL, MATSEQSELL, MATMPISELL 4371 M*/ 4372 4373 /*MC 4374 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4375 4376 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4377 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4378 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4379 for communicators controlling multiple processes. It is recommended that you call both of 4380 the above preallocation routines for simplicity. 4381 4382 Options Database Keys: 4383 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4384 4385 Level: beginner 4386 4387 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 4388 M*/ 4389 4390 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4391 #if defined(PETSC_HAVE_ELEMENTAL) 4392 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4393 #endif 4394 #if defined(PETSC_HAVE_SCALAPACK) 4395 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*); 4396 #endif 4397 #if defined(PETSC_HAVE_HYPRE) 4398 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4399 #endif 4400 4401 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4402 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4403 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4404 4405 /*@C 4406 MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored 4407 4408 Not Collective 4409 4410 Input Parameter: 4411 . mat - a MATSEQAIJ matrix 4412 4413 Output Parameter: 4414 . array - pointer to the data 4415 4416 Level: intermediate 4417 4418 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4419 @*/ 4420 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4421 { 4422 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4423 4424 PetscFunctionBegin; 4425 if (aij->ops->getarray) { 4426 PetscCall((*aij->ops->getarray)(A,array)); 4427 } else { 4428 *array = aij->a; 4429 } 4430 PetscFunctionReturn(0); 4431 } 4432 4433 /*@C 4434 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4435 4436 Not Collective 4437 4438 Input Parameters: 4439 + mat - a MATSEQAIJ matrix 4440 - array - pointer to the data 4441 4442 Level: intermediate 4443 4444 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4445 @*/ 4446 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4447 { 4448 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4449 4450 PetscFunctionBegin; 4451 if (aij->ops->restorearray) { 4452 PetscCall((*aij->ops->restorearray)(A,array)); 4453 } else { 4454 *array = NULL; 4455 } 4456 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4457 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4458 PetscFunctionReturn(0); 4459 } 4460 4461 /*@C 4462 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored 4463 4464 Not Collective 4465 4466 Input Parameter: 4467 . mat - a MATSEQAIJ matrix 4468 4469 Output Parameter: 4470 . array - pointer to the data 4471 4472 Level: intermediate 4473 4474 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead() 4475 @*/ 4476 PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array) 4477 { 4478 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4479 4480 PetscFunctionBegin; 4481 if (aij->ops->getarrayread) { 4482 PetscCall((*aij->ops->getarrayread)(A,array)); 4483 } else { 4484 *array = aij->a; 4485 } 4486 PetscFunctionReturn(0); 4487 } 4488 4489 /*@C 4490 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4491 4492 Not Collective 4493 4494 Input Parameter: 4495 . mat - a MATSEQAIJ matrix 4496 4497 Output Parameter: 4498 . array - pointer to the data 4499 4500 Level: intermediate 4501 4502 .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead() 4503 @*/ 4504 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array) 4505 { 4506 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4507 4508 PetscFunctionBegin; 4509 if (aij->ops->restorearrayread) { 4510 PetscCall((*aij->ops->restorearrayread)(A,array)); 4511 } else { 4512 *array = NULL; 4513 } 4514 PetscFunctionReturn(0); 4515 } 4516 4517 /*@C 4518 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a MATSEQAIJ matrix is stored 4519 4520 Not Collective 4521 4522 Input Parameter: 4523 . mat - a MATSEQAIJ matrix 4524 4525 Output Parameter: 4526 . array - pointer to the data 4527 4528 Level: intermediate 4529 4530 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead() 4531 @*/ 4532 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A,PetscScalar **array) 4533 { 4534 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4535 4536 PetscFunctionBegin; 4537 if (aij->ops->getarraywrite) { 4538 PetscCall((*aij->ops->getarraywrite)(A,array)); 4539 } else { 4540 *array = aij->a; 4541 } 4542 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4543 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4544 PetscFunctionReturn(0); 4545 } 4546 4547 /*@C 4548 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4549 4550 Not Collective 4551 4552 Input Parameter: 4553 . mat - a MATSEQAIJ matrix 4554 4555 Output Parameter: 4556 . array - pointer to the data 4557 4558 Level: intermediate 4559 4560 .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead() 4561 @*/ 4562 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A,PetscScalar **array) 4563 { 4564 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4565 4566 PetscFunctionBegin; 4567 if (aij->ops->restorearraywrite) { 4568 PetscCall((*aij->ops->restorearraywrite)(A,array)); 4569 } else { 4570 *array = NULL; 4571 } 4572 PetscFunctionReturn(0); 4573 } 4574 4575 /*@C 4576 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the SEQAIJ matrix 4577 4578 Not Collective 4579 4580 Input Parameter: 4581 . mat - a matrix of type MATSEQAIJ or its subclasses 4582 4583 Output Parameters: 4584 + i - row map array of the matrix 4585 . j - column index array of the matrix 4586 . a - data array of the matrix 4587 - memtype - memory type of the arrays 4588 4589 Notes: 4590 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4591 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4592 4593 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4594 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4595 4596 Level: Developer 4597 4598 .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead() 4599 @*/ 4600 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat,const PetscInt **i,const PetscInt **j,PetscScalar **a,PetscMemType *mtype) 4601 { 4602 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 4603 4604 PetscFunctionBegin; 4605 PetscCheck(mat->preallocated,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"matrix is not preallocated"); 4606 if (aij->ops->getcsrandmemtype) { 4607 PetscCall((*aij->ops->getcsrandmemtype)(mat,i,j,a,mtype)); 4608 } else { 4609 if (i) *i = aij->i; 4610 if (j) *j = aij->j; 4611 if (a) *a = aij->a; 4612 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4613 } 4614 PetscFunctionReturn(0); 4615 } 4616 4617 /*@C 4618 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4619 4620 Not Collective 4621 4622 Input Parameter: 4623 . mat - a MATSEQAIJ matrix 4624 4625 Output Parameter: 4626 . nz - the maximum number of nonzeros in any row 4627 4628 Level: intermediate 4629 4630 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4631 @*/ 4632 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4633 { 4634 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4635 4636 PetscFunctionBegin; 4637 *nz = aij->rmax; 4638 PetscFunctionReturn(0); 4639 } 4640 4641 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, const PetscInt coo_i[], const PetscInt coo_j[]) 4642 { 4643 MPI_Comm comm; 4644 PetscInt *i,*j; 4645 PetscInt M,N,row; 4646 PetscCount k,p,q,nneg,nnz,start,end; /* Index the coo array, so use PetscCount as their type */ 4647 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4648 PetscInt *Aj; 4649 PetscScalar *Aa; 4650 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ*)(mat->data); 4651 MatType rtype; 4652 PetscCount *perm,*jmap; 4653 4654 PetscFunctionBegin; 4655 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4656 PetscCall(PetscObjectGetComm((PetscObject)mat,&comm)); 4657 PetscCall(MatGetSize(mat,&M,&N)); 4658 PetscCall(PetscMalloc2(coo_n,&i,coo_n,&j)); 4659 PetscCall(PetscArraycpy(i,coo_i,coo_n)); /* Make a copy since we'll modify it */ 4660 PetscCall(PetscArraycpy(j,coo_j,coo_n)); 4661 PetscCall(PetscMalloc1(coo_n,&perm)); 4662 for (k=0; k<coo_n; k++) { /* Ignore entries with negative row or col indices */ 4663 if (j[k] < 0) i[k] = -1; 4664 perm[k] = k; 4665 } 4666 4667 /* Sort by row */ 4668 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n,i,j,perm)); 4669 for (k=0; k<coo_n; k++) {if (i[k] >= 0) break;} /* Advance k to the first row with a non-negative index */ 4670 nneg = k; 4671 PetscCall(PetscMalloc1(coo_n-nneg+1,&jmap)); /* +1 to make a CSR-like data structure. jmap[i] originally is the number of repeats for i-th nonzero */ 4672 nnz = 0; /* Total number of unique nonzeros to be counted */ 4673 jmap++; /* Inc jmap by 1 for convinience */ 4674 4675 PetscCall(PetscCalloc1(M+1,&Ai)); /* CSR of A */ 4676 PetscCall(PetscMalloc1(coo_n-nneg,&Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4677 4678 /* In each row, sort by column, then unique column indices to get row length */ 4679 Ai++; /* Inc by 1 for convinience */ 4680 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4681 while (k<coo_n) { 4682 row = i[k]; 4683 start = k; /* [start,end) indices for this row */ 4684 while (k<coo_n && i[k] == row) k++; 4685 end = k; 4686 PetscCall(PetscSortIntWithCountArray(end-start,j+start,perm+start)); 4687 /* Find number of unique col entries in this row */ 4688 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4689 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4690 Ai[row] = 1; 4691 nnz++; 4692 4693 for (p=start+1; p<end; p++) { /* Scan remaining nonzero in this row */ 4694 if (j[p] != j[p-1]) { /* Meet a new nonzero */ 4695 q++; 4696 jmap[q] = 1; 4697 Aj[q] = j[p]; 4698 Ai[row]++; 4699 nnz++; 4700 } else { 4701 jmap[q]++; 4702 } 4703 } 4704 q++; /* Move to next row and thus next unique nonzero */ 4705 } 4706 PetscCall(PetscFree2(i,j)); 4707 4708 Ai--; /* Back to the beginning of Ai[] */ 4709 for (k=0; k<M; k++) Ai[k+1] += Ai[k]; 4710 jmap--; /* Back to the beginning of jmap[] */ 4711 jmap[0] = 0; 4712 for (k=0; k<nnz; k++) jmap[k+1] += jmap[k]; 4713 if (nnz < coo_n-nneg) { /* Realloc with actual number of unique nonzeros */ 4714 PetscCount *jmap_new; 4715 PetscInt *Aj_new; 4716 4717 PetscCall(PetscMalloc1(nnz+1,&jmap_new)); 4718 PetscCall(PetscArraycpy(jmap_new,jmap,nnz+1)); 4719 PetscCall(PetscFree(jmap)); 4720 jmap = jmap_new; 4721 4722 PetscCall(PetscMalloc1(nnz,&Aj_new)); 4723 PetscCall(PetscArraycpy(Aj_new,Aj,nnz)); 4724 PetscCall(PetscFree(Aj)); 4725 Aj = Aj_new; 4726 } 4727 4728 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4729 PetscCount *perm_new; 4730 4731 PetscCall(PetscMalloc1(coo_n-nneg,&perm_new)); 4732 PetscCall(PetscArraycpy(perm_new,perm+nneg,coo_n-nneg)); 4733 PetscCall(PetscFree(perm)); 4734 perm = perm_new; 4735 } 4736 4737 PetscCall(MatGetRootType_Private(mat,&rtype)); 4738 PetscCall(PetscCalloc1(nnz,&Aa)); /* Zero the matrix */ 4739 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF,M,N,Ai,Aj,Aa,rtype,mat)); 4740 4741 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4742 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4743 /* Record COO fields */ 4744 seqaij->coo_n = coo_n; 4745 seqaij->Atot = coo_n-nneg; /* Annz is seqaij->nz, so no need to record that again */ 4746 seqaij->jmap = jmap; /* of length nnz+1 */ 4747 seqaij->perm = perm; 4748 PetscFunctionReturn(0); 4749 } 4750 4751 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A,const PetscScalar v[],InsertMode imode) 4752 { 4753 Mat_SeqAIJ *aseq = (Mat_SeqAIJ*)A->data; 4754 PetscCount i,j,Annz = aseq->nz; 4755 PetscCount *perm = aseq->perm,*jmap = aseq->jmap; 4756 PetscScalar *Aa; 4757 4758 PetscFunctionBegin; 4759 PetscCall(MatSeqAIJGetArray(A,&Aa)); 4760 for (i=0; i<Annz; i++) { 4761 PetscScalar sum = 0.0; 4762 for (j=jmap[i]; j<jmap[i+1]; j++) sum += v[perm[j]]; 4763 Aa[i] = (imode == INSERT_VALUES? 0.0 : Aa[i]) + sum; 4764 } 4765 PetscCall(MatSeqAIJRestoreArray(A,&Aa)); 4766 PetscFunctionReturn(0); 4767 } 4768 4769 #if defined(PETSC_HAVE_CUDA) 4770 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat,MatType,MatReuse,Mat*); 4771 #endif 4772 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4773 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat,MatType,MatReuse,Mat*); 4774 #endif 4775 4776 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4777 { 4778 Mat_SeqAIJ *b; 4779 PetscMPIInt size; 4780 4781 PetscFunctionBegin; 4782 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 4783 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4784 4785 PetscCall(PetscNewLog(B,&b)); 4786 4787 B->data = (void*)b; 4788 4789 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 4790 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4791 4792 b->row = NULL; 4793 b->col = NULL; 4794 b->icol = NULL; 4795 b->reallocs = 0; 4796 b->ignorezeroentries = PETSC_FALSE; 4797 b->roworiented = PETSC_TRUE; 4798 b->nonew = 0; 4799 b->diag = NULL; 4800 b->solve_work = NULL; 4801 B->spptr = NULL; 4802 b->saved_values = NULL; 4803 b->idiag = NULL; 4804 b->mdiag = NULL; 4805 b->ssor_work = NULL; 4806 b->omega = 1.0; 4807 b->fshift = 0.0; 4808 b->idiagvalid = PETSC_FALSE; 4809 b->ibdiagvalid = PETSC_FALSE; 4810 b->keepnonzeropattern = PETSC_FALSE; 4811 4812 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4813 4814 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4815 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ)); 4816 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ)); 4817 #endif 4818 4819 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ)); 4820 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ)); 4821 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ)); 4822 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ)); 4823 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ)); 4824 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM)); 4825 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL)); 4826 #if defined(PETSC_HAVE_MKL_SPARSE) 4827 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL)); 4828 #endif 4829 #if defined(PETSC_HAVE_CUDA) 4830 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4831 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4832 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",MatProductSetFromOptions_SeqAIJ)); 4833 #endif 4834 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4835 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijkokkos_C",MatConvert_SeqAIJ_SeqAIJKokkos)); 4836 #endif 4837 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL)); 4838 #if defined(PETSC_HAVE_ELEMENTAL) 4839 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental)); 4840 #endif 4841 #if defined(PETSC_HAVE_SCALAPACK) 4842 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK)); 4843 #endif 4844 #if defined(PETSC_HAVE_HYPRE) 4845 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE)); 4846 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4847 #endif 4848 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense)); 4849 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL)); 4850 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS)); 4851 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ)); 4852 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ)); 4853 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ)); 4854 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ)); 4855 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4856 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ)); 4857 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ)); 4858 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ)); 4859 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4860 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJKron_C",MatSeqAIJKron_SeqAIJ)); 4861 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_SeqAIJ)); 4862 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_SeqAIJ)); 4863 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4864 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4865 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4866 PetscFunctionReturn(0); 4867 } 4868 4869 /* 4870 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4871 */ 4872 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4873 { 4874 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data,*a = (Mat_SeqAIJ*)A->data; 4875 PetscInt m = A->rmap->n,i; 4876 4877 PetscFunctionBegin; 4878 PetscCheckFalse(!A->assembled && cpvalues!=MAT_DO_NOT_COPY_VALUES,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot duplicate unassembled matrix"); 4879 4880 C->factortype = A->factortype; 4881 c->row = NULL; 4882 c->col = NULL; 4883 c->icol = NULL; 4884 c->reallocs = 0; 4885 4886 C->assembled = A->assembled; 4887 C->preallocated = A->preallocated; 4888 4889 if (A->preallocated) { 4890 PetscCall(PetscLayoutReference(A->rmap,&C->rmap)); 4891 PetscCall(PetscLayoutReference(A->cmap,&C->cmap)); 4892 4893 PetscCall(PetscMalloc1(m,&c->imax)); 4894 PetscCall(PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt))); 4895 PetscCall(PetscMalloc1(m,&c->ilen)); 4896 PetscCall(PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt))); 4897 PetscCall(PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt))); 4898 4899 /* allocate the matrix space */ 4900 if (mallocmatspace) { 4901 PetscCall(PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i)); 4902 PetscCall(PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt))); 4903 4904 c->singlemalloc = PETSC_TRUE; 4905 4906 PetscCall(PetscArraycpy(c->i,a->i,m+1)); 4907 if (m > 0) { 4908 PetscCall(PetscArraycpy(c->j,a->j,a->i[m])); 4909 if (cpvalues == MAT_COPY_VALUES) { 4910 const PetscScalar *aa; 4911 4912 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4913 PetscCall(PetscArraycpy(c->a,aa,a->i[m])); 4914 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4915 } else { 4916 PetscCall(PetscArrayzero(c->a,a->i[m])); 4917 } 4918 } 4919 } 4920 4921 c->ignorezeroentries = a->ignorezeroentries; 4922 c->roworiented = a->roworiented; 4923 c->nonew = a->nonew; 4924 if (a->diag) { 4925 PetscCall(PetscMalloc1(m+1,&c->diag)); 4926 PetscCall(PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt))); 4927 PetscCall(PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt))); 4928 } else c->diag = NULL; 4929 4930 c->solve_work = NULL; 4931 c->saved_values = NULL; 4932 c->idiag = NULL; 4933 c->ssor_work = NULL; 4934 c->keepnonzeropattern = a->keepnonzeropattern; 4935 c->free_a = PETSC_TRUE; 4936 c->free_ij = PETSC_TRUE; 4937 4938 c->rmax = a->rmax; 4939 c->nz = a->nz; 4940 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4941 4942 c->compressedrow.use = a->compressedrow.use; 4943 c->compressedrow.nrows = a->compressedrow.nrows; 4944 if (a->compressedrow.use) { 4945 i = a->compressedrow.nrows; 4946 PetscCall(PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex)); 4947 PetscCall(PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1)); 4948 PetscCall(PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i)); 4949 } else { 4950 c->compressedrow.use = PETSC_FALSE; 4951 c->compressedrow.i = NULL; 4952 c->compressedrow.rindex = NULL; 4953 } 4954 c->nonzerorowcnt = a->nonzerorowcnt; 4955 C->nonzerostate = A->nonzerostate; 4956 4957 PetscCall(MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C)); 4958 } 4959 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist)); 4960 PetscFunctionReturn(0); 4961 } 4962 4963 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4964 { 4965 PetscFunctionBegin; 4966 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 4967 PetscCall(MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 4968 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4969 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 4970 } 4971 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 4972 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE)); 4973 PetscFunctionReturn(0); 4974 } 4975 4976 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4977 { 4978 PetscBool isbinary, ishdf5; 4979 4980 PetscFunctionBegin; 4981 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4982 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4983 /* force binary viewer to load .info file if it has not yet done so */ 4984 PetscCall(PetscViewerSetUp(viewer)); 4985 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 4986 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5)); 4987 if (isbinary) { 4988 PetscCall(MatLoad_SeqAIJ_Binary(newMat,viewer)); 4989 } else if (ishdf5) { 4990 #if defined(PETSC_HAVE_HDF5) 4991 PetscCall(MatLoad_AIJ_HDF5(newMat,viewer)); 4992 #else 4993 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4994 #endif 4995 } else { 4996 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 4997 } 4998 PetscFunctionReturn(0); 4999 } 5000 5001 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 5002 { 5003 Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data; 5004 PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i; 5005 5006 PetscFunctionBegin; 5007 PetscCall(PetscViewerSetUp(viewer)); 5008 5009 /* read in matrix header */ 5010 PetscCall(PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT)); 5011 PetscCheck(header[0] == MAT_FILE_CLASSID,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 5012 M = header[1]; N = header[2]; nz = header[3]; 5013 PetscCheck(M >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%" PetscInt_FMT ") in file is negative",M); 5014 PetscCheck(N >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%" PetscInt_FMT ") in file is negative",N); 5015 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 5016 5017 /* set block sizes from the viewer's .info file */ 5018 PetscCall(MatLoad_Binary_BlockSizes(mat,viewer)); 5019 /* set local and global sizes if not set already */ 5020 if (mat->rmap->n < 0) mat->rmap->n = M; 5021 if (mat->cmap->n < 0) mat->cmap->n = N; 5022 if (mat->rmap->N < 0) mat->rmap->N = M; 5023 if (mat->cmap->N < 0) mat->cmap->N = N; 5024 PetscCall(PetscLayoutSetUp(mat->rmap)); 5025 PetscCall(PetscLayoutSetUp(mat->cmap)); 5026 5027 /* check if the matrix sizes are correct */ 5028 PetscCall(MatGetSize(mat,&rows,&cols)); 5029 PetscCheckFalse(M != rows || N != cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5030 5031 /* read in row lengths */ 5032 PetscCall(PetscMalloc1(M,&rowlens)); 5033 PetscCall(PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT)); 5034 /* check if sum(rowlens) is same as nz */ 5035 sum = 0; for (i=0; i<M; i++) sum += rowlens[i]; 5036 PetscCheck(sum == nz,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT,nz,sum); 5037 /* preallocate and check sizes */ 5038 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens)); 5039 PetscCall(MatGetSize(mat,&rows,&cols)); 5040 PetscCheckFalse(M != rows || N != cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5041 /* store row lengths */ 5042 PetscCall(PetscArraycpy(a->ilen,rowlens,M)); 5043 PetscCall(PetscFree(rowlens)); 5044 5045 /* fill in "i" row pointers */ 5046 a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i]; 5047 /* read in "j" column indices */ 5048 PetscCall(PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT)); 5049 /* read in "a" nonzero values */ 5050 PetscCall(PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR)); 5051 5052 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 5053 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 5054 PetscFunctionReturn(0); 5055 } 5056 5057 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 5058 { 5059 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 5060 const PetscScalar *aa,*ba; 5061 #if defined(PETSC_USE_COMPLEX) 5062 PetscInt k; 5063 #endif 5064 5065 PetscFunctionBegin; 5066 /* If the matrix dimensions are not equal,or no of nonzeros */ 5067 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 5068 *flg = PETSC_FALSE; 5069 PetscFunctionReturn(0); 5070 } 5071 5072 /* if the a->i are the same */ 5073 PetscCall(PetscArraycmp(a->i,b->i,A->rmap->n+1,flg)); 5074 if (!*flg) PetscFunctionReturn(0); 5075 5076 /* if a->j are the same */ 5077 PetscCall(PetscArraycmp(a->j,b->j,a->nz,flg)); 5078 if (!*flg) PetscFunctionReturn(0); 5079 5080 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 5081 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 5082 /* if a->a are the same */ 5083 #if defined(PETSC_USE_COMPLEX) 5084 for (k=0; k<a->nz; k++) { 5085 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5086 *flg = PETSC_FALSE; 5087 PetscFunctionReturn(0); 5088 } 5089 } 5090 #else 5091 PetscCall(PetscArraycmp(aa,ba,a->nz,flg)); 5092 #endif 5093 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 5094 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 5095 PetscFunctionReturn(0); 5096 } 5097 5098 /*@ 5099 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 5100 provided by the user. 5101 5102 Collective 5103 5104 Input Parameters: 5105 + comm - must be an MPI communicator of size 1 5106 . m - number of rows 5107 . n - number of columns 5108 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5109 . j - column indices 5110 - a - matrix values 5111 5112 Output Parameter: 5113 . mat - the matrix 5114 5115 Level: intermediate 5116 5117 Notes: 5118 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5119 once the matrix is destroyed and not before 5120 5121 You cannot set new nonzero locations into this matrix, that will generate an error. 5122 5123 The i and j indices are 0 based 5124 5125 The format which is used for the sparse matrix input, is equivalent to a 5126 row-major ordering.. i.e for the following matrix, the input data expected is 5127 as shown 5128 5129 $ 1 0 0 5130 $ 2 0 3 5131 $ 4 5 6 5132 $ 5133 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5134 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5135 $ v = {1,2,3,4,5,6} [size = 6] 5136 5137 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 5138 5139 @*/ 5140 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 5141 { 5142 PetscInt ii; 5143 Mat_SeqAIJ *aij; 5144 PetscInt jj; 5145 5146 PetscFunctionBegin; 5147 PetscCheckFalse(m > 0 && i[0],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 5148 PetscCall(MatCreate(comm,mat)); 5149 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5150 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5151 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5152 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL)); 5153 aij = (Mat_SeqAIJ*)(*mat)->data; 5154 PetscCall(PetscMalloc1(m,&aij->imax)); 5155 PetscCall(PetscMalloc1(m,&aij->ilen)); 5156 5157 aij->i = i; 5158 aij->j = j; 5159 aij->a = a; 5160 aij->singlemalloc = PETSC_FALSE; 5161 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5162 aij->free_a = PETSC_FALSE; 5163 aij->free_ij = PETSC_FALSE; 5164 5165 for (ii=0,aij->nonzerorowcnt=0,aij->rmax=0; ii<m; ii++) { 5166 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 5167 if (PetscDefined(USE_DEBUG)) { 5168 PetscCheckFalse(i[ii+1] - i[ii] < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT,ii,i[ii+1] - i[ii]); 5169 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 5170 PetscCheck(j[jj] >= j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is not sorted",jj-i[ii],j[jj],ii); 5171 PetscCheck(j[jj] != j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is identical to previous entry",jj-i[ii],j[jj],ii); 5172 } 5173 } 5174 } 5175 if (PetscDefined(USE_DEBUG)) { 5176 for (ii=0; ii<aij->i[m]; ii++) { 5177 PetscCheck(j[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5178 PetscCheckFalse(j[ii] > n - 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5179 } 5180 } 5181 5182 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5183 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5184 PetscFunctionReturn(0); 5185 } 5186 5187 /*@C 5188 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 5189 provided by the user. 5190 5191 Collective 5192 5193 Input Parameters: 5194 + comm - must be an MPI communicator of size 1 5195 . m - number of rows 5196 . n - number of columns 5197 . i - row indices 5198 . j - column indices 5199 . a - matrix values 5200 . nz - number of nonzeros 5201 - idx - 0 or 1 based 5202 5203 Output Parameter: 5204 . mat - the matrix 5205 5206 Level: intermediate 5207 5208 Notes: 5209 The i and j indices are 0 based. The format which is used for the sparse matrix input, is equivalent to a row-major ordering. i.e for the following matrix, 5210 the input data expected is as shown 5211 .vb 5212 1 0 0 5213 2 0 3 5214 4 5 6 5215 5216 i = {0,1,1,2,2,2} 5217 j = {0,0,2,0,1,2} 5218 v = {1,2,3,4,5,6} 5219 .ve 5220 5221 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 5222 5223 @*/ 5224 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 5225 { 5226 PetscInt ii, *nnz, one = 1,row,col; 5227 5228 PetscFunctionBegin; 5229 PetscCall(PetscCalloc1(m,&nnz)); 5230 for (ii = 0; ii < nz; ii++) { 5231 nnz[i[ii] - !!idx] += 1; 5232 } 5233 PetscCall(MatCreate(comm,mat)); 5234 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5235 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5236 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz)); 5237 for (ii = 0; ii < nz; ii++) { 5238 if (idx) { 5239 row = i[ii] - 1; 5240 col = j[ii] - 1; 5241 } else { 5242 row = i[ii]; 5243 col = j[ii]; 5244 } 5245 PetscCall(MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES)); 5246 } 5247 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5248 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5249 PetscCall(PetscFree(nnz)); 5250 PetscFunctionReturn(0); 5251 } 5252 5253 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5254 { 5255 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 5256 5257 PetscFunctionBegin; 5258 a->idiagvalid = PETSC_FALSE; 5259 a->ibdiagvalid = PETSC_FALSE; 5260 5261 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5262 PetscFunctionReturn(0); 5263 } 5264 5265 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 5266 { 5267 PetscMPIInt size; 5268 5269 PetscFunctionBegin; 5270 PetscCallMPI(MPI_Comm_size(comm,&size)); 5271 if (size == 1) { 5272 if (scall == MAT_INITIAL_MATRIX) { 5273 PetscCall(MatDuplicate(inmat,MAT_COPY_VALUES,outmat)); 5274 } else { 5275 PetscCall(MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN)); 5276 } 5277 } else { 5278 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat)); 5279 } 5280 PetscFunctionReturn(0); 5281 } 5282 5283 /* 5284 Permute A into C's *local* index space using rowemb,colemb. 5285 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5286 of [0,m), colemb is in [0,n). 5287 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5288 */ 5289 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 5290 { 5291 /* If making this function public, change the error returned in this function away from _PLIB. */ 5292 Mat_SeqAIJ *Baij; 5293 PetscBool seqaij; 5294 PetscInt m,n,*nz,i,j,count; 5295 PetscScalar v; 5296 const PetscInt *rowindices,*colindices; 5297 5298 PetscFunctionBegin; 5299 if (!B) PetscFunctionReturn(0); 5300 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5301 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij)); 5302 PetscCheck(seqaij,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 5303 if (rowemb) { 5304 PetscCall(ISGetLocalSize(rowemb,&m)); 5305 PetscCheck(m == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Row IS of size %" PetscInt_FMT " is incompatible with matrix row size %" PetscInt_FMT,m,B->rmap->n); 5306 } else { 5307 PetscCheck(C->rmap->n == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 5308 } 5309 if (colemb) { 5310 PetscCall(ISGetLocalSize(colemb,&n)); 5311 PetscCheck(n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diag col IS of size %" PetscInt_FMT " is incompatible with input matrix col size %" PetscInt_FMT,n,B->cmap->n); 5312 } else { 5313 PetscCheck(C->cmap->n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 5314 } 5315 5316 Baij = (Mat_SeqAIJ*)(B->data); 5317 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5318 PetscCall(PetscMalloc1(B->rmap->n,&nz)); 5319 for (i=0; i<B->rmap->n; i++) { 5320 nz[i] = Baij->i[i+1] - Baij->i[i]; 5321 } 5322 PetscCall(MatSeqAIJSetPreallocation(C,0,nz)); 5323 PetscCall(PetscFree(nz)); 5324 } 5325 if (pattern == SUBSET_NONZERO_PATTERN) { 5326 PetscCall(MatZeroEntries(C)); 5327 } 5328 count = 0; 5329 rowindices = NULL; 5330 colindices = NULL; 5331 if (rowemb) { 5332 PetscCall(ISGetIndices(rowemb,&rowindices)); 5333 } 5334 if (colemb) { 5335 PetscCall(ISGetIndices(colemb,&colindices)); 5336 } 5337 for (i=0; i<B->rmap->n; i++) { 5338 PetscInt row; 5339 row = i; 5340 if (rowindices) row = rowindices[i]; 5341 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 5342 PetscInt col; 5343 col = Baij->j[count]; 5344 if (colindices) col = colindices[col]; 5345 v = Baij->a[count]; 5346 PetscCall(MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES)); 5347 ++count; 5348 } 5349 } 5350 /* FIXME: set C's nonzerostate correctly. */ 5351 /* Assembly for C is necessary. */ 5352 C->preallocated = PETSC_TRUE; 5353 C->assembled = PETSC_TRUE; 5354 C->was_assembled = PETSC_FALSE; 5355 PetscFunctionReturn(0); 5356 } 5357 5358 PetscFunctionList MatSeqAIJList = NULL; 5359 5360 /*@C 5361 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 5362 5363 Collective on Mat 5364 5365 Input Parameters: 5366 + mat - the matrix object 5367 - matype - matrix type 5368 5369 Options Database Key: 5370 . -mat_seqai_type <method> - for example seqaijcrl 5371 5372 Level: intermediate 5373 5374 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 5375 @*/ 5376 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5377 { 5378 PetscBool sametype; 5379 PetscErrorCode (*r)(Mat,MatType,MatReuse,Mat*); 5380 5381 PetscFunctionBegin; 5382 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 5383 PetscCall(PetscObjectTypeCompare((PetscObject)mat,matype,&sametype)); 5384 if (sametype) PetscFunctionReturn(0); 5385 5386 PetscCall(PetscFunctionListFind(MatSeqAIJList,matype,&r)); 5387 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 5388 PetscCall((*r)(mat,matype,MAT_INPLACE_MATRIX,&mat)); 5389 PetscFunctionReturn(0); 5390 } 5391 5392 /*@C 5393 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 5394 5395 Not Collective 5396 5397 Input Parameters: 5398 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 5399 - function - routine to convert to subtype 5400 5401 Notes: 5402 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 5403 5404 Then, your matrix can be chosen with the procedural interface at runtime via the option 5405 $ -mat_seqaij_type my_mat 5406 5407 Level: advanced 5408 5409 .seealso: MatSeqAIJRegisterAll() 5410 5411 Level: advanced 5412 @*/ 5413 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 5414 { 5415 PetscFunctionBegin; 5416 PetscCall(MatInitializePackage()); 5417 PetscCall(PetscFunctionListAdd(&MatSeqAIJList,sname,function)); 5418 PetscFunctionReturn(0); 5419 } 5420 5421 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5422 5423 /*@C 5424 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 5425 5426 Not Collective 5427 5428 Level: advanced 5429 5430 .seealso: MatRegisterAll(), MatSeqAIJRegister() 5431 @*/ 5432 PetscErrorCode MatSeqAIJRegisterAll(void) 5433 { 5434 PetscFunctionBegin; 5435 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5436 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5437 5438 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5439 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5440 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5441 #if defined(PETSC_HAVE_MKL_SPARSE) 5442 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5443 #endif 5444 #if defined(PETSC_HAVE_CUDA) 5445 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5446 #endif 5447 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5448 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5449 #endif 5450 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5451 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5452 #endif 5453 PetscFunctionReturn(0); 5454 } 5455 5456 /* 5457 Special version for direct calls from Fortran 5458 */ 5459 #include <petsc/private/fortranimpl.h> 5460 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5461 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5462 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5463 #define matsetvaluesseqaij_ matsetvaluesseqaij 5464 #endif 5465 5466 /* Change these macros so can be used in void function */ 5467 5468 /* Change these macros so can be used in void function */ 5469 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5470 #undef PetscCall 5471 #define PetscCall(...) do { \ 5472 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5473 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5474 *_ierr = PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr_msv_mpiaij,PETSC_ERROR_REPEAT," "); \ 5475 return; \ 5476 } \ 5477 } while (0) 5478 5479 #undef SETERRQ 5480 #define SETERRQ(comm,ierr,...) do { \ 5481 *_ierr = PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__); \ 5482 return; \ 5483 } while (0) 5484 5485 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 5486 { 5487 Mat A = *AA; 5488 PetscInt m = *mm, n = *nn; 5489 InsertMode is = *isis; 5490 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 5491 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 5492 PetscInt *imax,*ai,*ailen; 5493 PetscInt *aj,nonew = a->nonew,lastcol = -1; 5494 MatScalar *ap,value,*aa; 5495 PetscBool ignorezeroentries = a->ignorezeroentries; 5496 PetscBool roworiented = a->roworiented; 5497 5498 PetscFunctionBegin; 5499 MatCheckPreallocated(A,1); 5500 imax = a->imax; 5501 ai = a->i; 5502 ailen = a->ilen; 5503 aj = a->j; 5504 aa = a->a; 5505 5506 for (k=0; k<m; k++) { /* loop over added rows */ 5507 row = im[k]; 5508 if (row < 0) continue; 5509 PetscCheck(row < A->rmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 5510 rp = aj + ai[row]; ap = aa + ai[row]; 5511 rmax = imax[row]; nrow = ailen[row]; 5512 low = 0; 5513 high = nrow; 5514 for (l=0; l<n; l++) { /* loop over added columns */ 5515 if (in[l] < 0) continue; 5516 PetscCheck(in[l] < A->cmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 5517 col = in[l]; 5518 if (roworiented) value = v[l + k*n]; 5519 else value = v[k + l*m]; 5520 5521 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5522 5523 if (col <= lastcol) low = 0; 5524 else high = nrow; 5525 lastcol = col; 5526 while (high-low > 5) { 5527 t = (low+high)/2; 5528 if (rp[t] > col) high = t; 5529 else low = t; 5530 } 5531 for (i=low; i<high; i++) { 5532 if (rp[i] > col) break; 5533 if (rp[i] == col) { 5534 if (is == ADD_VALUES) ap[i] += value; 5535 else ap[i] = value; 5536 goto noinsert; 5537 } 5538 } 5539 if (value == 0.0 && ignorezeroentries) goto noinsert; 5540 if (nonew == 1) goto noinsert; 5541 PetscCheck(nonew != -1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 5542 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 5543 N = nrow++ - 1; a->nz++; high++; 5544 /* shift up all the later entries in this row */ 5545 for (ii=N; ii>=i; ii--) { 5546 rp[ii+1] = rp[ii]; 5547 ap[ii+1] = ap[ii]; 5548 } 5549 rp[i] = col; 5550 ap[i] = value; 5551 A->nonzerostate++; 5552 noinsert:; 5553 low = i + 1; 5554 } 5555 ailen[row] = nrow; 5556 } 5557 PetscFunctionReturnVoid(); 5558 } 5559 /* Undefining these here since they were redefined from their original definition above! No 5560 * other PETSc functions should be defined past this point, as it is impossible to recover the 5561 * original definitions */ 5562 #undef PetscCall 5563 #undef SETERRQ 5564