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