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