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