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