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) { 469 ap[i] += value; 470 (void)PetscLogFlops(1.0); 471 } 472 else ap[i] = value; 473 } 474 low = i + 1; 475 goto noinsert; 476 } 477 } 478 if (value == 0.0 && ignorezeroentries && row != col) goto noinsert; 479 if (nonew == 1) goto noinsert; 480 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 481 if (A->structure_only) { 482 MatSeqXAIJReallocateAIJ_structure_only(A,A->rmap->n,1,nrow,row,col,rmax,ai,aj,rp,imax,nonew,MatScalar); 483 } else { 484 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 485 } 486 N = nrow++ - 1; a->nz++; high++; 487 /* shift up all the later entries in this row */ 488 for (ii=N; ii>=i; ii--) { 489 rp[ii+1] = rp[ii]; 490 if (!A->structure_only) ap[ii+1] = ap[ii]; 491 } 492 rp[i] = col; 493 if (!A->structure_only) ap[i] = value; 494 low = i + 1; 495 A->nonzerostate++; 496 noinsert:; 497 } 498 ailen[row] = nrow; 499 } 500 PetscFunctionReturn(0); 501 } 502 503 504 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 505 { 506 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 507 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 508 PetscInt *ai = a->i,*ailen = a->ilen; 509 MatScalar *ap,*aa = a->a; 510 511 PetscFunctionBegin; 512 for (k=0; k<m; k++) { /* loop over rows */ 513 row = im[k]; 514 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 515 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); 516 rp = aj + ai[row]; ap = aa + ai[row]; 517 nrow = ailen[row]; 518 for (l=0; l<n; l++) { /* loop over columns */ 519 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 520 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); 521 col = in[l]; 522 high = nrow; low = 0; /* assume unsorted */ 523 while (high-low > 5) { 524 t = (low+high)/2; 525 if (rp[t] > col) high = t; 526 else low = t; 527 } 528 for (i=low; i<high; i++) { 529 if (rp[i] > col) break; 530 if (rp[i] == col) { 531 *v++ = ap[i]; 532 goto finished; 533 } 534 } 535 *v++ = 0.0; 536 finished:; 537 } 538 } 539 PetscFunctionReturn(0); 540 } 541 542 543 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 544 { 545 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 546 PetscErrorCode ierr; 547 PetscInt i,*col_lens; 548 int fd; 549 FILE *file; 550 551 PetscFunctionBegin; 552 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 553 ierr = PetscMalloc1(4+A->rmap->n,&col_lens);CHKERRQ(ierr); 554 555 col_lens[0] = MAT_FILE_CLASSID; 556 col_lens[1] = A->rmap->n; 557 col_lens[2] = A->cmap->n; 558 col_lens[3] = a->nz; 559 560 /* store lengths of each row and write (including header) to file */ 561 for (i=0; i<A->rmap->n; i++) { 562 col_lens[4+i] = a->i[i+1] - a->i[i]; 563 } 564 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 565 ierr = PetscFree(col_lens);CHKERRQ(ierr); 566 567 /* store column indices (zero start index) */ 568 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 569 570 /* store nonzero values */ 571 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 572 573 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 574 if (file) { 575 fprintf(file,"-matload_block_size %d\n",(int)PetscAbs(A->rmap->bs)); 576 } 577 PetscFunctionReturn(0); 578 } 579 580 static PetscErrorCode MatView_SeqAIJ_ASCII_structonly(Mat A,PetscViewer viewer) 581 { 582 PetscErrorCode ierr; 583 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 584 PetscInt i,k,m=A->rmap->N; 585 586 PetscFunctionBegin; 587 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 588 for (i=0; i<m; i++) { 589 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 590 for (k=a->i[i]; k<a->i[i+1]; k++) { 591 ierr = PetscViewerASCIIPrintf(viewer," (%D) ",a->j[k]);CHKERRQ(ierr); 592 } 593 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 594 } 595 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 596 PetscFunctionReturn(0); 597 } 598 599 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 600 601 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 602 { 603 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 604 PetscErrorCode ierr; 605 PetscInt i,j,m = A->rmap->n; 606 const char *name; 607 PetscViewerFormat format; 608 609 PetscFunctionBegin; 610 if (A->structure_only) { 611 ierr = MatView_SeqAIJ_ASCII_structonly(A,viewer);CHKERRQ(ierr); 612 PetscFunctionReturn(0); 613 } 614 615 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 616 if (format == PETSC_VIEWER_ASCII_MATLAB) { 617 PetscInt nofinalvalue = 0; 618 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) { 619 /* Need a dummy value to ensure the dimension of the matrix. */ 620 nofinalvalue = 1; 621 } 622 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 623 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 624 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 625 #if defined(PETSC_USE_COMPLEX) 626 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,4);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 627 #else 628 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 629 #endif 630 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 631 632 for (i=0; i<m; i++) { 633 for (j=a->i[i]; j<a->i[i+1]; j++) { 634 #if defined(PETSC_USE_COMPLEX) 635 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); 636 #else 637 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+1,(double)a->a[j]);CHKERRQ(ierr); 638 #endif 639 } 640 } 641 if (nofinalvalue) { 642 #if defined(PETSC_USE_COMPLEX) 643 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e %18.16e\n",m,A->cmap->n,0.,0.);CHKERRQ(ierr); 644 #else 645 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 646 #endif 647 } 648 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 649 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 650 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 651 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 652 PetscFunctionReturn(0); 653 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 654 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 655 for (i=0; i<m; i++) { 656 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 657 for (j=a->i[i]; j<a->i[i+1]; j++) { 658 #if defined(PETSC_USE_COMPLEX) 659 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 660 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 661 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 662 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 663 } else if (PetscRealPart(a->a[j]) != 0.0) { 664 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 665 } 666 #else 667 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr);} 668 #endif 669 } 670 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 671 } 672 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 673 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 674 PetscInt nzd=0,fshift=1,*sptr; 675 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 676 ierr = PetscMalloc1(m+1,&sptr);CHKERRQ(ierr); 677 for (i=0; i<m; i++) { 678 sptr[i] = nzd+1; 679 for (j=a->i[i]; j<a->i[i+1]; j++) { 680 if (a->j[j] >= i) { 681 #if defined(PETSC_USE_COMPLEX) 682 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 683 #else 684 if (a->a[j] != 0.0) nzd++; 685 #endif 686 } 687 } 688 } 689 sptr[m] = nzd+1; 690 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 691 for (i=0; i<m+1; i+=6) { 692 if (i+4<m) { 693 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); 694 } else if (i+3<m) { 695 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); 696 } else if (i+2<m) { 697 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr); 698 } else if (i+1<m) { 699 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr); 700 } else if (i<m) { 701 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr); 702 } else { 703 ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr); 704 } 705 } 706 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 707 ierr = PetscFree(sptr);CHKERRQ(ierr); 708 for (i=0; i<m; i++) { 709 for (j=a->i[i]; j<a->i[i+1]; j++) { 710 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 711 } 712 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 713 } 714 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 715 for (i=0; i<m; i++) { 716 for (j=a->i[i]; j<a->i[i+1]; j++) { 717 if (a->j[j] >= i) { 718 #if defined(PETSC_USE_COMPLEX) 719 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 720 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 721 } 722 #else 723 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j]);CHKERRQ(ierr);} 724 #endif 725 } 726 } 727 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 728 } 729 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 730 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 731 PetscInt cnt = 0,jcnt; 732 PetscScalar value; 733 #if defined(PETSC_USE_COMPLEX) 734 PetscBool realonly = PETSC_TRUE; 735 736 for (i=0; i<a->i[m]; i++) { 737 if (PetscImaginaryPart(a->a[i]) != 0.0) { 738 realonly = PETSC_FALSE; 739 break; 740 } 741 } 742 #endif 743 744 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 745 for (i=0; i<m; i++) { 746 jcnt = 0; 747 for (j=0; j<A->cmap->n; j++) { 748 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 749 value = a->a[cnt++]; 750 jcnt++; 751 } else { 752 value = 0.0; 753 } 754 #if defined(PETSC_USE_COMPLEX) 755 if (realonly) { 756 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value));CHKERRQ(ierr); 757 } else { 758 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value));CHKERRQ(ierr); 759 } 760 #else 761 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value);CHKERRQ(ierr); 762 #endif 763 } 764 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 765 } 766 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 767 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 768 PetscInt fshift=1; 769 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 770 #if defined(PETSC_USE_COMPLEX) 771 ierr = PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n");CHKERRQ(ierr); 772 #else 773 ierr = PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n");CHKERRQ(ierr); 774 #endif 775 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 776 for (i=0; i<m; i++) { 777 for (j=a->i[i]; j<a->i[i+1]; j++) { 778 #if defined(PETSC_USE_COMPLEX) 779 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); 780 #else 781 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j]);CHKERRQ(ierr); 782 #endif 783 } 784 } 785 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 786 } else { 787 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 788 if (A->factortype) { 789 for (i=0; i<m; i++) { 790 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 791 /* L part */ 792 for (j=a->i[i]; j<a->i[i+1]; j++) { 793 #if defined(PETSC_USE_COMPLEX) 794 if (PetscImaginaryPart(a->a[j]) > 0.0) { 795 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 796 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 797 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));CHKERRQ(ierr); 798 } else { 799 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 800 } 801 #else 802 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 803 #endif 804 } 805 /* diagonal */ 806 j = a->diag[i]; 807 #if defined(PETSC_USE_COMPLEX) 808 if (PetscImaginaryPart(a->a[j]) > 0.0) { 809 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); 810 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 811 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); 812 } else { 813 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]));CHKERRQ(ierr); 814 } 815 #else 816 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)(1.0/a->a[j]));CHKERRQ(ierr); 817 #endif 818 819 /* U part */ 820 for (j=a->diag[i+1]+1; j<a->diag[i]; j++) { 821 #if defined(PETSC_USE_COMPLEX) 822 if (PetscImaginaryPart(a->a[j]) > 0.0) { 823 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 824 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 825 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));CHKERRQ(ierr); 826 } else { 827 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 828 } 829 #else 830 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 831 #endif 832 } 833 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 834 } 835 } else { 836 for (i=0; i<m; i++) { 837 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 838 for (j=a->i[i]; j<a->i[i+1]; j++) { 839 #if defined(PETSC_USE_COMPLEX) 840 if (PetscImaginaryPart(a->a[j]) > 0.0) { 841 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 842 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 843 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 844 } else { 845 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 846 } 847 #else 848 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 849 #endif 850 } 851 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 852 } 853 } 854 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 855 } 856 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 857 PetscFunctionReturn(0); 858 } 859 860 #include <petscdraw.h> 861 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 862 { 863 Mat A = (Mat) Aa; 864 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 865 PetscErrorCode ierr; 866 PetscInt i,j,m = A->rmap->n; 867 int color; 868 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 869 PetscViewer viewer; 870 PetscViewerFormat format; 871 872 PetscFunctionBegin; 873 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 874 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 875 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 876 877 /* loop over matrix elements drawing boxes */ 878 879 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 880 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 881 /* Blue for negative, Cyan for zero and Red for positive */ 882 color = PETSC_DRAW_BLUE; 883 for (i=0; i<m; i++) { 884 y_l = m - i - 1.0; y_r = y_l + 1.0; 885 for (j=a->i[i]; j<a->i[i+1]; j++) { 886 x_l = a->j[j]; x_r = x_l + 1.0; 887 if (PetscRealPart(a->a[j]) >= 0.) continue; 888 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 889 } 890 } 891 color = PETSC_DRAW_CYAN; 892 for (i=0; i<m; i++) { 893 y_l = m - i - 1.0; y_r = y_l + 1.0; 894 for (j=a->i[i]; j<a->i[i+1]; j++) { 895 x_l = a->j[j]; x_r = x_l + 1.0; 896 if (a->a[j] != 0.) continue; 897 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 898 } 899 } 900 color = PETSC_DRAW_RED; 901 for (i=0; i<m; i++) { 902 y_l = m - i - 1.0; y_r = y_l + 1.0; 903 for (j=a->i[i]; j<a->i[i+1]; j++) { 904 x_l = a->j[j]; x_r = x_l + 1.0; 905 if (PetscRealPart(a->a[j]) <= 0.) continue; 906 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 907 } 908 } 909 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 910 } else { 911 /* use contour shading to indicate magnitude of values */ 912 /* first determine max of all nonzero values */ 913 PetscReal minv = 0.0, maxv = 0.0; 914 PetscInt nz = a->nz, count = 0; 915 PetscDraw popup; 916 917 for (i=0; i<nz; i++) { 918 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 919 } 920 if (minv >= maxv) maxv = minv + PETSC_SMALL; 921 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 922 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 923 924 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 925 for (i=0; i<m; i++) { 926 y_l = m - i - 1.0; 927 y_r = y_l + 1.0; 928 for (j=a->i[i]; j<a->i[i+1]; j++) { 929 x_l = a->j[j]; 930 x_r = x_l + 1.0; 931 color = PetscDrawRealToColor(PetscAbsScalar(a->a[count]),minv,maxv); 932 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 933 count++; 934 } 935 } 936 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 937 } 938 PetscFunctionReturn(0); 939 } 940 941 #include <petscdraw.h> 942 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 943 { 944 PetscErrorCode ierr; 945 PetscDraw draw; 946 PetscReal xr,yr,xl,yl,h,w; 947 PetscBool isnull; 948 949 PetscFunctionBegin; 950 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 951 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 952 if (isnull) PetscFunctionReturn(0); 953 954 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 955 xr += w; yr += h; xl = -w; yl = -h; 956 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 957 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 958 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 959 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 960 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 961 PetscFunctionReturn(0); 962 } 963 964 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 965 { 966 PetscErrorCode ierr; 967 PetscBool iascii,isbinary,isdraw; 968 969 PetscFunctionBegin; 970 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 971 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 972 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 973 if (iascii) { 974 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 975 } else if (isbinary) { 976 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 977 } else if (isdraw) { 978 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 979 } 980 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 981 PetscFunctionReturn(0); 982 } 983 984 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 985 { 986 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 987 PetscErrorCode ierr; 988 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 989 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 990 MatScalar *aa = a->a,*ap; 991 PetscReal ratio = 0.6; 992 993 PetscFunctionBegin; 994 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 995 996 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 997 for (i=1; i<m; i++) { 998 /* move each row back by the amount of empty slots (fshift) before it*/ 999 fshift += imax[i-1] - ailen[i-1]; 1000 rmax = PetscMax(rmax,ailen[i]); 1001 if (fshift) { 1002 ip = aj + ai[i]; 1003 ap = aa + ai[i]; 1004 N = ailen[i]; 1005 for (j=0; j<N; j++) { 1006 ip[j-fshift] = ip[j]; 1007 if (!A->structure_only) ap[j-fshift] = ap[j]; 1008 } 1009 } 1010 ai[i] = ai[i-1] + ailen[i-1]; 1011 } 1012 if (m) { 1013 fshift += imax[m-1] - ailen[m-1]; 1014 ai[m] = ai[m-1] + ailen[m-1]; 1015 } 1016 1017 /* reset ilen and imax for each row */ 1018 a->nonzerorowcnt = 0; 1019 if (A->structure_only) { 1020 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 1021 } else { /* !A->structure_only */ 1022 for (i=0; i<m; i++) { 1023 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1024 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 1025 } 1026 } 1027 a->nz = ai[m]; 1028 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); 1029 1030 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1031 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 1032 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 1033 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 1034 1035 A->info.mallocs += a->reallocs; 1036 a->reallocs = 0; 1037 A->info.nz_unneeded = (PetscReal)fshift; 1038 a->rmax = rmax; 1039 1040 if (!A->structure_only) { 1041 ierr = MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 1042 } 1043 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 1044 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1045 PetscFunctionReturn(0); 1046 } 1047 1048 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1049 { 1050 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1051 PetscInt i,nz = a->nz; 1052 MatScalar *aa = a->a; 1053 PetscErrorCode ierr; 1054 1055 PetscFunctionBegin; 1056 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1057 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1058 PetscFunctionReturn(0); 1059 } 1060 1061 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1062 { 1063 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1064 PetscInt i,nz = a->nz; 1065 MatScalar *aa = a->a; 1066 PetscErrorCode ierr; 1067 1068 PetscFunctionBegin; 1069 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1070 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1071 PetscFunctionReturn(0); 1072 } 1073 1074 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1075 { 1076 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1077 PetscErrorCode ierr; 1078 1079 PetscFunctionBegin; 1080 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 1081 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1082 PetscFunctionReturn(0); 1083 } 1084 1085 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1086 { 1087 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1088 PetscErrorCode ierr; 1089 1090 PetscFunctionBegin; 1091 #if defined(PETSC_USE_LOG) 1092 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 1093 #endif 1094 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 1095 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1096 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1097 ierr = PetscFree(a->diag);CHKERRQ(ierr); 1098 ierr = PetscFree(a->ibdiag);CHKERRQ(ierr); 1099 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 1100 ierr = PetscFree(a->ipre);CHKERRQ(ierr); 1101 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 1102 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 1103 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 1104 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 1105 ierr = ISColoringDestroy(&a->coloring);CHKERRQ(ierr); 1106 ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);CHKERRQ(ierr); 1107 ierr = PetscFree(a->matmult_abdense);CHKERRQ(ierr); 1108 1109 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 1110 ierr = PetscFree(A->data);CHKERRQ(ierr); 1111 1112 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1113 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL);CHKERRQ(ierr); 1114 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);CHKERRQ(ierr); 1115 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);CHKERRQ(ierr); 1116 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL);CHKERRQ(ierr); 1117 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL);CHKERRQ(ierr); 1118 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL);CHKERRQ(ierr); 1119 #if defined(PETSC_HAVE_ELEMENTAL) 1120 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL);CHKERRQ(ierr); 1121 #endif 1122 #if defined(PETSC_HAVE_HYPRE) 1123 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_hypre_C",NULL);CHKERRQ(ierr); 1124 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMatMatMult_transpose_seqaij_seqaij_C",NULL);CHKERRQ(ierr); 1125 #endif 1126 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1127 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsell_C",NULL);CHKERRQ(ierr); 1128 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_is_C",NULL);CHKERRQ(ierr); 1129 ierr = PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);CHKERRQ(ierr); 1130 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);CHKERRQ(ierr); 1131 ierr = PetscObjectComposeFunction((PetscObject)A,"MatResetPreallocation_C",NULL);CHKERRQ(ierr); 1132 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr); 1133 ierr = PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);CHKERRQ(ierr); 1134 ierr = PetscObjectComposeFunction((PetscObject)A,"MatPtAP_is_seqaij_C",NULL);CHKERRQ(ierr); 1135 PetscFunctionReturn(0); 1136 } 1137 1138 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1139 { 1140 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1141 PetscErrorCode ierr; 1142 1143 PetscFunctionBegin; 1144 switch (op) { 1145 case MAT_ROW_ORIENTED: 1146 a->roworiented = flg; 1147 break; 1148 case MAT_KEEP_NONZERO_PATTERN: 1149 a->keepnonzeropattern = flg; 1150 break; 1151 case MAT_NEW_NONZERO_LOCATIONS: 1152 a->nonew = (flg ? 0 : 1); 1153 break; 1154 case MAT_NEW_NONZERO_LOCATION_ERR: 1155 a->nonew = (flg ? -1 : 0); 1156 break; 1157 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1158 a->nonew = (flg ? -2 : 0); 1159 break; 1160 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1161 a->nounused = (flg ? -1 : 0); 1162 break; 1163 case MAT_IGNORE_ZERO_ENTRIES: 1164 a->ignorezeroentries = flg; 1165 break; 1166 case MAT_SPD: 1167 case MAT_SYMMETRIC: 1168 case MAT_STRUCTURALLY_SYMMETRIC: 1169 case MAT_HERMITIAN: 1170 case MAT_SYMMETRY_ETERNAL: 1171 case MAT_STRUCTURE_ONLY: 1172 /* These options are handled directly by MatSetOption() */ 1173 break; 1174 case MAT_NEW_DIAGONALS: 1175 case MAT_IGNORE_OFF_PROC_ENTRIES: 1176 case MAT_USE_HASH_TABLE: 1177 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1178 break; 1179 case MAT_USE_INODES: 1180 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 1181 break; 1182 case MAT_SUBMAT_SINGLEIS: 1183 A->submat_singleis = flg; 1184 break; 1185 default: 1186 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1187 } 1188 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 1189 PetscFunctionReturn(0); 1190 } 1191 1192 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1193 { 1194 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1195 PetscErrorCode ierr; 1196 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 1197 PetscScalar *aa=a->a,*x,zero=0.0; 1198 1199 PetscFunctionBegin; 1200 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1201 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1202 1203 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1204 PetscInt *diag=a->diag; 1205 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1206 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1207 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1208 PetscFunctionReturn(0); 1209 } 1210 1211 ierr = VecSet(v,zero);CHKERRQ(ierr); 1212 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1213 for (i=0; i<n; i++) { 1214 nz = ai[i+1] - ai[i]; 1215 if (!nz) x[i] = 0.0; 1216 for (j=ai[i]; j<ai[i+1]; j++) { 1217 if (aj[j] == i) { 1218 x[i] = aa[j]; 1219 break; 1220 } 1221 } 1222 } 1223 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1224 PetscFunctionReturn(0); 1225 } 1226 1227 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1228 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1229 { 1230 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1231 PetscScalar *y; 1232 const PetscScalar *x; 1233 PetscErrorCode ierr; 1234 PetscInt m = A->rmap->n; 1235 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1236 const MatScalar *v; 1237 PetscScalar alpha; 1238 PetscInt n,i,j; 1239 const PetscInt *idx,*ii,*ridx=NULL; 1240 Mat_CompressedRow cprow = a->compressedrow; 1241 PetscBool usecprow = cprow.use; 1242 #endif 1243 1244 PetscFunctionBegin; 1245 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 1246 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1247 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1248 1249 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1250 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 1251 #else 1252 if (usecprow) { 1253 m = cprow.nrows; 1254 ii = cprow.i; 1255 ridx = cprow.rindex; 1256 } else { 1257 ii = a->i; 1258 } 1259 for (i=0; i<m; i++) { 1260 idx = a->j + ii[i]; 1261 v = a->a + ii[i]; 1262 n = ii[i+1] - ii[i]; 1263 if (usecprow) { 1264 alpha = x[ridx[i]]; 1265 } else { 1266 alpha = x[i]; 1267 } 1268 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1269 } 1270 #endif 1271 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1272 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1273 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1274 PetscFunctionReturn(0); 1275 } 1276 1277 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1278 { 1279 PetscErrorCode ierr; 1280 1281 PetscFunctionBegin; 1282 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 1283 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 1284 PetscFunctionReturn(0); 1285 } 1286 1287 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1288 1289 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1290 { 1291 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1292 PetscScalar *y; 1293 const PetscScalar *x; 1294 const MatScalar *aa; 1295 PetscErrorCode ierr; 1296 PetscInt m=A->rmap->n; 1297 const PetscInt *aj,*ii,*ridx=NULL; 1298 PetscInt n,i; 1299 PetscScalar sum; 1300 PetscBool usecprow=a->compressedrow.use; 1301 1302 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1303 #pragma disjoint(*x,*y,*aa) 1304 #endif 1305 1306 PetscFunctionBegin; 1307 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1308 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1309 ii = a->i; 1310 if (usecprow) { /* use compressed row format */ 1311 ierr = PetscMemzero(y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1312 m = a->compressedrow.nrows; 1313 ii = a->compressedrow.i; 1314 ridx = a->compressedrow.rindex; 1315 for (i=0; i<m; i++) { 1316 n = ii[i+1] - ii[i]; 1317 aj = a->j + ii[i]; 1318 aa = a->a + ii[i]; 1319 sum = 0.0; 1320 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1321 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1322 y[*ridx++] = sum; 1323 } 1324 } else { /* do not use compressed row format */ 1325 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1326 aj = a->j; 1327 aa = a->a; 1328 fortranmultaij_(&m,x,ii,aj,aa,y); 1329 #else 1330 for (i=0; i<m; i++) { 1331 n = ii[i+1] - ii[i]; 1332 aj = a->j + ii[i]; 1333 aa = a->a + ii[i]; 1334 sum = 0.0; 1335 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1336 y[i] = sum; 1337 } 1338 #endif 1339 } 1340 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1341 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1342 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1343 PetscFunctionReturn(0); 1344 } 1345 1346 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1347 { 1348 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1349 PetscScalar *y; 1350 const PetscScalar *x; 1351 const MatScalar *aa; 1352 PetscErrorCode ierr; 1353 PetscInt m=A->rmap->n; 1354 const PetscInt *aj,*ii,*ridx=NULL; 1355 PetscInt n,i,nonzerorow=0; 1356 PetscScalar sum; 1357 PetscBool usecprow=a->compressedrow.use; 1358 1359 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1360 #pragma disjoint(*x,*y,*aa) 1361 #endif 1362 1363 PetscFunctionBegin; 1364 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1365 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1366 if (usecprow) { /* use compressed row format */ 1367 m = a->compressedrow.nrows; 1368 ii = a->compressedrow.i; 1369 ridx = a->compressedrow.rindex; 1370 for (i=0; i<m; i++) { 1371 n = ii[i+1] - ii[i]; 1372 aj = a->j + ii[i]; 1373 aa = a->a + ii[i]; 1374 sum = 0.0; 1375 nonzerorow += (n>0); 1376 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1377 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1378 y[*ridx++] = sum; 1379 } 1380 } else { /* do not use compressed row format */ 1381 ii = a->i; 1382 for (i=0; i<m; i++) { 1383 n = ii[i+1] - ii[i]; 1384 aj = a->j + ii[i]; 1385 aa = a->a + ii[i]; 1386 sum = 0.0; 1387 nonzerorow += (n>0); 1388 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1389 y[i] = sum; 1390 } 1391 } 1392 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1393 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1394 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1395 PetscFunctionReturn(0); 1396 } 1397 1398 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1399 { 1400 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1401 PetscScalar *y,*z; 1402 const PetscScalar *x; 1403 const MatScalar *aa; 1404 PetscErrorCode ierr; 1405 PetscInt m = A->rmap->n,*aj,*ii; 1406 PetscInt n,i,*ridx=NULL; 1407 PetscScalar sum; 1408 PetscBool usecprow=a->compressedrow.use; 1409 1410 PetscFunctionBegin; 1411 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1412 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1413 if (usecprow) { /* use compressed row format */ 1414 if (zz != yy) { 1415 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1416 } 1417 m = a->compressedrow.nrows; 1418 ii = a->compressedrow.i; 1419 ridx = a->compressedrow.rindex; 1420 for (i=0; i<m; i++) { 1421 n = ii[i+1] - ii[i]; 1422 aj = a->j + ii[i]; 1423 aa = a->a + ii[i]; 1424 sum = y[*ridx]; 1425 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1426 z[*ridx++] = sum; 1427 } 1428 } else { /* do not use compressed row format */ 1429 ii = a->i; 1430 for (i=0; i<m; i++) { 1431 n = ii[i+1] - ii[i]; 1432 aj = a->j + ii[i]; 1433 aa = a->a + ii[i]; 1434 sum = y[i]; 1435 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1436 z[i] = sum; 1437 } 1438 } 1439 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1440 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1441 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1442 PetscFunctionReturn(0); 1443 } 1444 1445 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1446 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1447 { 1448 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1449 PetscScalar *y,*z; 1450 const PetscScalar *x; 1451 const MatScalar *aa; 1452 PetscErrorCode ierr; 1453 const PetscInt *aj,*ii,*ridx=NULL; 1454 PetscInt m = A->rmap->n,n,i; 1455 PetscScalar sum; 1456 PetscBool usecprow=a->compressedrow.use; 1457 1458 PetscFunctionBegin; 1459 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1460 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1461 if (usecprow) { /* use compressed row format */ 1462 if (zz != yy) { 1463 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1464 } 1465 m = a->compressedrow.nrows; 1466 ii = a->compressedrow.i; 1467 ridx = a->compressedrow.rindex; 1468 for (i=0; i<m; i++) { 1469 n = ii[i+1] - ii[i]; 1470 aj = a->j + ii[i]; 1471 aa = a->a + ii[i]; 1472 sum = y[*ridx]; 1473 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1474 z[*ridx++] = sum; 1475 } 1476 } else { /* do not use compressed row format */ 1477 ii = a->i; 1478 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1479 aj = a->j; 1480 aa = a->a; 1481 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1482 #else 1483 for (i=0; i<m; i++) { 1484 n = ii[i+1] - ii[i]; 1485 aj = a->j + ii[i]; 1486 aa = a->a + ii[i]; 1487 sum = y[i]; 1488 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1489 z[i] = sum; 1490 } 1491 #endif 1492 } 1493 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1494 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1495 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1496 PetscFunctionReturn(0); 1497 } 1498 1499 /* 1500 Adds diagonal pointers to sparse matrix structure. 1501 */ 1502 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1503 { 1504 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1505 PetscErrorCode ierr; 1506 PetscInt i,j,m = A->rmap->n; 1507 1508 PetscFunctionBegin; 1509 if (!a->diag) { 1510 ierr = PetscMalloc1(m,&a->diag);CHKERRQ(ierr); 1511 ierr = PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));CHKERRQ(ierr); 1512 } 1513 for (i=0; i<A->rmap->n; i++) { 1514 a->diag[i] = a->i[i+1]; 1515 for (j=a->i[i]; j<a->i[i+1]; j++) { 1516 if (a->j[j] == i) { 1517 a->diag[i] = j; 1518 break; 1519 } 1520 } 1521 } 1522 PetscFunctionReturn(0); 1523 } 1524 1525 PetscErrorCode MatShift_SeqAIJ(Mat A,PetscScalar v) 1526 { 1527 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1528 const PetscInt *diag = (const PetscInt*)a->diag; 1529 const PetscInt *ii = (const PetscInt*) a->i; 1530 PetscInt i,*mdiag = NULL; 1531 PetscErrorCode ierr; 1532 PetscInt cnt = 0; /* how many diagonals are missing */ 1533 1534 PetscFunctionBegin; 1535 if (!A->preallocated || !a->nz) { 1536 ierr = MatSeqAIJSetPreallocation(A,1,NULL);CHKERRQ(ierr); 1537 ierr = MatShift_Basic(A,v);CHKERRQ(ierr); 1538 PetscFunctionReturn(0); 1539 } 1540 1541 if (a->diagonaldense) { 1542 cnt = 0; 1543 } else { 1544 ierr = PetscCalloc1(A->rmap->n,&mdiag);CHKERRQ(ierr); 1545 for (i=0; i<A->rmap->n; i++) { 1546 if (diag[i] >= ii[i+1]) { 1547 cnt++; 1548 mdiag[i] = 1; 1549 } 1550 } 1551 } 1552 if (!cnt) { 1553 ierr = MatShift_Basic(A,v);CHKERRQ(ierr); 1554 } else { 1555 PetscScalar *olda = a->a; /* preserve pointers to current matrix nonzeros structure and values */ 1556 PetscInt *oldj = a->j, *oldi = a->i; 1557 PetscBool singlemalloc = a->singlemalloc,free_a = a->free_a,free_ij = a->free_ij; 1558 1559 a->a = NULL; 1560 a->j = NULL; 1561 a->i = NULL; 1562 /* increase the values in imax for each row where a diagonal is being inserted then reallocate the matrix data structures */ 1563 for (i=0; i<A->rmap->n; i++) { 1564 a->imax[i] += mdiag[i]; 1565 a->imax[i] = PetscMin(a->imax[i],A->cmap->n); 1566 } 1567 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax);CHKERRQ(ierr); 1568 1569 /* copy old values into new matrix data structure */ 1570 for (i=0; i<A->rmap->n; i++) { 1571 ierr = MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES);CHKERRQ(ierr); 1572 if (i < A->cmap->n) { 1573 ierr = MatSetValue(A,i,i,v,ADD_VALUES);CHKERRQ(ierr); 1574 } 1575 } 1576 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1577 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1578 if (singlemalloc) { 1579 ierr = PetscFree3(olda,oldj,oldi);CHKERRQ(ierr); 1580 } else { 1581 if (free_a) {ierr = PetscFree(olda);CHKERRQ(ierr);} 1582 if (free_ij) {ierr = PetscFree(oldj);CHKERRQ(ierr);} 1583 if (free_ij) {ierr = PetscFree(oldi);CHKERRQ(ierr);} 1584 } 1585 } 1586 ierr = PetscFree(mdiag);CHKERRQ(ierr); 1587 a->diagonaldense = PETSC_TRUE; 1588 PetscFunctionReturn(0); 1589 } 1590 1591 /* 1592 Checks for missing diagonals 1593 */ 1594 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1595 { 1596 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1597 PetscInt *diag,*ii = a->i,i; 1598 PetscErrorCode ierr; 1599 1600 PetscFunctionBegin; 1601 *missing = PETSC_FALSE; 1602 if (A->rmap->n > 0 && !ii) { 1603 *missing = PETSC_TRUE; 1604 if (d) *d = 0; 1605 ierr = PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");CHKERRQ(ierr); 1606 } else { 1607 diag = a->diag; 1608 for (i=0; i<A->rmap->n; i++) { 1609 if (diag[i] >= ii[i+1]) { 1610 *missing = PETSC_TRUE; 1611 if (d) *d = i; 1612 ierr = PetscInfo1(A,"Matrix is missing diagonal number %D\n",i);CHKERRQ(ierr); 1613 break; 1614 } 1615 } 1616 } 1617 PetscFunctionReturn(0); 1618 } 1619 1620 #include <petscblaslapack.h> 1621 #include <petsc/private/kernels/blockinvert.h> 1622 1623 /* 1624 Note that values is allocated externally by the PC and then passed into this routine 1625 */ 1626 PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag) 1627 { 1628 PetscErrorCode ierr; 1629 PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots; 1630 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 1631 const PetscReal shift = 0.0; 1632 PetscInt ipvt[5]; 1633 PetscScalar work[25],*v_work; 1634 1635 PetscFunctionBegin; 1636 allowzeropivot = PetscNot(A->erroriffailure); 1637 for (i=0; i<nblocks; i++) ncnt += bsizes[i]; 1638 if (ncnt != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Total blocksizes %D doesn't match number matrix rows %D",ncnt,n); 1639 for (i=0; i<nblocks; i++) { 1640 bsizemax = PetscMax(bsizemax,bsizes[i]); 1641 } 1642 ierr = PetscMalloc1(bsizemax,&indx);CHKERRQ(ierr); 1643 if (bsizemax > 7) { 1644 ierr = PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots);CHKERRQ(ierr); 1645 } 1646 ncnt = 0; 1647 for (i=0; i<nblocks; i++) { 1648 for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j; 1649 ierr = MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag);CHKERRQ(ierr); 1650 switch (bsizes[i]) { 1651 case 1: 1652 *diag = 1.0/(*diag); 1653 break; 1654 case 2: 1655 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1656 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1657 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 1658 break; 1659 case 3: 1660 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1661 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1662 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 1663 break; 1664 case 4: 1665 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1666 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1667 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 1668 break; 1669 case 5: 1670 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1671 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1672 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 1673 break; 1674 case 6: 1675 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1676 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1677 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 1678 break; 1679 case 7: 1680 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1681 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1682 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 1683 break; 1684 default: 1685 ierr = PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1686 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1687 ierr = PetscKernel_A_gets_transpose_A_N(diag,bsizes[i]);CHKERRQ(ierr); 1688 } 1689 ncnt += bsizes[i]; 1690 diag += bsizes[i]*bsizes[i]; 1691 } 1692 if (bsizemax > 7) { 1693 ierr = PetscFree2(v_work,v_pivots);CHKERRQ(ierr); 1694 } 1695 ierr = PetscFree(indx);CHKERRQ(ierr); 1696 PetscFunctionReturn(0); 1697 } 1698 1699 /* 1700 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1701 */ 1702 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1703 { 1704 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1705 PetscErrorCode ierr; 1706 PetscInt i,*diag,m = A->rmap->n; 1707 MatScalar *v = a->a; 1708 PetscScalar *idiag,*mdiag; 1709 1710 PetscFunctionBegin; 1711 if (a->idiagvalid) PetscFunctionReturn(0); 1712 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1713 diag = a->diag; 1714 if (!a->idiag) { 1715 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1716 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1717 v = a->a; 1718 } 1719 mdiag = a->mdiag; 1720 idiag = a->idiag; 1721 1722 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1723 for (i=0; i<m; i++) { 1724 mdiag[i] = v[diag[i]]; 1725 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1726 if (PetscRealPart(fshift)) { 1727 ierr = PetscInfo1(A,"Zero diagonal on row %D\n",i);CHKERRQ(ierr); 1728 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1729 A->factorerror_zeropivot_value = 0.0; 1730 A->factorerror_zeropivot_row = i; 1731 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1732 } 1733 idiag[i] = 1.0/v[diag[i]]; 1734 } 1735 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1736 } else { 1737 for (i=0; i<m; i++) { 1738 mdiag[i] = v[diag[i]]; 1739 idiag[i] = omega/(fshift + v[diag[i]]); 1740 } 1741 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1742 } 1743 a->idiagvalid = PETSC_TRUE; 1744 PetscFunctionReturn(0); 1745 } 1746 1747 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1748 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1749 { 1750 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1751 PetscScalar *x,d,sum,*t,scale; 1752 const MatScalar *v,*idiag=0,*mdiag; 1753 const PetscScalar *b, *bs,*xb, *ts; 1754 PetscErrorCode ierr; 1755 PetscInt n,m = A->rmap->n,i; 1756 const PetscInt *idx,*diag; 1757 1758 PetscFunctionBegin; 1759 its = its*lits; 1760 1761 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1762 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1763 a->fshift = fshift; 1764 a->omega = omega; 1765 1766 diag = a->diag; 1767 t = a->ssor_work; 1768 idiag = a->idiag; 1769 mdiag = a->mdiag; 1770 1771 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1772 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1773 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1774 if (flag == SOR_APPLY_UPPER) { 1775 /* apply (U + D/omega) to the vector */ 1776 bs = b; 1777 for (i=0; i<m; i++) { 1778 d = fshift + mdiag[i]; 1779 n = a->i[i+1] - diag[i] - 1; 1780 idx = a->j + diag[i] + 1; 1781 v = a->a + diag[i] + 1; 1782 sum = b[i]*d/omega; 1783 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1784 x[i] = sum; 1785 } 1786 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1787 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1788 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1789 PetscFunctionReturn(0); 1790 } 1791 1792 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1793 else if (flag & SOR_EISENSTAT) { 1794 /* Let A = L + U + D; where L is lower trianglar, 1795 U is upper triangular, E = D/omega; This routine applies 1796 1797 (L + E)^{-1} A (U + E)^{-1} 1798 1799 to a vector efficiently using Eisenstat's trick. 1800 */ 1801 scale = (2.0/omega) - 1.0; 1802 1803 /* x = (E + U)^{-1} b */ 1804 for (i=m-1; i>=0; i--) { 1805 n = a->i[i+1] - diag[i] - 1; 1806 idx = a->j + diag[i] + 1; 1807 v = a->a + diag[i] + 1; 1808 sum = b[i]; 1809 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1810 x[i] = sum*idiag[i]; 1811 } 1812 1813 /* t = b - (2*E - D)x */ 1814 v = a->a; 1815 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1816 1817 /* t = (E + L)^{-1}t */ 1818 ts = t; 1819 diag = a->diag; 1820 for (i=0; i<m; i++) { 1821 n = diag[i] - a->i[i]; 1822 idx = a->j + a->i[i]; 1823 v = a->a + a->i[i]; 1824 sum = t[i]; 1825 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1826 t[i] = sum*idiag[i]; 1827 /* x = x + t */ 1828 x[i] += t[i]; 1829 } 1830 1831 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1832 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1833 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1834 PetscFunctionReturn(0); 1835 } 1836 if (flag & SOR_ZERO_INITIAL_GUESS) { 1837 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1838 for (i=0; i<m; i++) { 1839 n = diag[i] - a->i[i]; 1840 idx = a->j + a->i[i]; 1841 v = a->a + a->i[i]; 1842 sum = b[i]; 1843 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1844 t[i] = sum; 1845 x[i] = sum*idiag[i]; 1846 } 1847 xb = t; 1848 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1849 } else xb = b; 1850 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1851 for (i=m-1; i>=0; i--) { 1852 n = a->i[i+1] - diag[i] - 1; 1853 idx = a->j + diag[i] + 1; 1854 v = a->a + diag[i] + 1; 1855 sum = xb[i]; 1856 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1857 if (xb == b) { 1858 x[i] = sum*idiag[i]; 1859 } else { 1860 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1861 } 1862 } 1863 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1864 } 1865 its--; 1866 } 1867 while (its--) { 1868 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1869 for (i=0; i<m; i++) { 1870 /* lower */ 1871 n = diag[i] - a->i[i]; 1872 idx = a->j + a->i[i]; 1873 v = a->a + a->i[i]; 1874 sum = b[i]; 1875 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1876 t[i] = sum; /* save application of the lower-triangular part */ 1877 /* upper */ 1878 n = a->i[i+1] - diag[i] - 1; 1879 idx = a->j + diag[i] + 1; 1880 v = a->a + diag[i] + 1; 1881 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1882 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1883 } 1884 xb = t; 1885 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1886 } else xb = b; 1887 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1888 for (i=m-1; i>=0; i--) { 1889 sum = xb[i]; 1890 if (xb == b) { 1891 /* whole matrix (no checkpointing available) */ 1892 n = a->i[i+1] - a->i[i]; 1893 idx = a->j + a->i[i]; 1894 v = a->a + a->i[i]; 1895 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1896 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1897 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1898 n = a->i[i+1] - diag[i] - 1; 1899 idx = a->j + diag[i] + 1; 1900 v = a->a + diag[i] + 1; 1901 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1902 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1903 } 1904 } 1905 if (xb == b) { 1906 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1907 } else { 1908 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1909 } 1910 } 1911 } 1912 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1913 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1914 PetscFunctionReturn(0); 1915 } 1916 1917 1918 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1919 { 1920 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1921 1922 PetscFunctionBegin; 1923 info->block_size = 1.0; 1924 info->nz_allocated = (double)a->maxnz; 1925 info->nz_used = (double)a->nz; 1926 info->nz_unneeded = (double)(a->maxnz - a->nz); 1927 info->assemblies = (double)A->num_ass; 1928 info->mallocs = (double)A->info.mallocs; 1929 info->memory = ((PetscObject)A)->mem; 1930 if (A->factortype) { 1931 info->fill_ratio_given = A->info.fill_ratio_given; 1932 info->fill_ratio_needed = A->info.fill_ratio_needed; 1933 info->factor_mallocs = A->info.factor_mallocs; 1934 } else { 1935 info->fill_ratio_given = 0; 1936 info->fill_ratio_needed = 0; 1937 info->factor_mallocs = 0; 1938 } 1939 PetscFunctionReturn(0); 1940 } 1941 1942 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1943 { 1944 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1945 PetscInt i,m = A->rmap->n - 1; 1946 PetscErrorCode ierr; 1947 const PetscScalar *xx; 1948 PetscScalar *bb; 1949 PetscInt d = 0; 1950 1951 PetscFunctionBegin; 1952 if (x && b) { 1953 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1954 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1955 for (i=0; i<N; i++) { 1956 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1957 if (rows[i] >= A->cmap->n) continue; 1958 bb[rows[i]] = diag*xx[rows[i]]; 1959 } 1960 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1961 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1962 } 1963 1964 if (a->keepnonzeropattern) { 1965 for (i=0; i<N; i++) { 1966 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1967 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1968 } 1969 if (diag != 0.0) { 1970 for (i=0; i<N; i++) { 1971 d = rows[i]; 1972 if (rows[i] >= A->cmap->n) continue; 1973 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); 1974 } 1975 for (i=0; i<N; i++) { 1976 if (rows[i] >= A->cmap->n) continue; 1977 a->a[a->diag[rows[i]]] = diag; 1978 } 1979 } 1980 } else { 1981 if (diag != 0.0) { 1982 for (i=0; i<N; i++) { 1983 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1984 if (a->ilen[rows[i]] > 0) { 1985 if (rows[i] >= A->cmap->n) { 1986 a->ilen[rows[i]] = 0; 1987 } else { 1988 a->ilen[rows[i]] = 1; 1989 a->a[a->i[rows[i]]] = diag; 1990 a->j[a->i[rows[i]]] = rows[i]; 1991 } 1992 } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */ 1993 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1994 } 1995 } 1996 } else { 1997 for (i=0; i<N; i++) { 1998 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1999 a->ilen[rows[i]] = 0; 2000 } 2001 } 2002 A->nonzerostate++; 2003 } 2004 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2005 PetscFunctionReturn(0); 2006 } 2007 2008 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2009 { 2010 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2011 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2012 PetscErrorCode ierr; 2013 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2014 const PetscScalar *xx; 2015 PetscScalar *bb; 2016 2017 PetscFunctionBegin; 2018 if (x && b) { 2019 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2020 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 2021 vecs = PETSC_TRUE; 2022 } 2023 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 2024 for (i=0; i<N; i++) { 2025 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2026 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 2027 2028 zeroed[rows[i]] = PETSC_TRUE; 2029 } 2030 for (i=0; i<A->rmap->n; i++) { 2031 if (!zeroed[i]) { 2032 for (j=a->i[i]; j<a->i[i+1]; j++) { 2033 if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) { 2034 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 2035 a->a[j] = 0.0; 2036 } 2037 } 2038 } else if (vecs && i < A->cmap->N) bb[i] = diag*xx[i]; 2039 } 2040 if (x && b) { 2041 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2042 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 2043 } 2044 ierr = PetscFree(zeroed);CHKERRQ(ierr); 2045 if (diag != 0.0) { 2046 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 2047 if (missing) { 2048 for (i=0; i<N; i++) { 2049 if (rows[i] >= A->cmap->N) continue; 2050 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]); 2051 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 2052 } 2053 } else { 2054 for (i=0; i<N; i++) { 2055 a->a[a->diag[rows[i]]] = diag; 2056 } 2057 } 2058 } 2059 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2060 PetscFunctionReturn(0); 2061 } 2062 2063 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2064 { 2065 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2066 PetscInt *itmp; 2067 2068 PetscFunctionBegin; 2069 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 2070 2071 *nz = a->i[row+1] - a->i[row]; 2072 if (v) *v = a->a + a->i[row]; 2073 if (idx) { 2074 itmp = a->j + a->i[row]; 2075 if (*nz) *idx = itmp; 2076 else *idx = 0; 2077 } 2078 PetscFunctionReturn(0); 2079 } 2080 2081 /* remove this function? */ 2082 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2083 { 2084 PetscFunctionBegin; 2085 PetscFunctionReturn(0); 2086 } 2087 2088 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2089 { 2090 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2091 MatScalar *v = a->a; 2092 PetscReal sum = 0.0; 2093 PetscErrorCode ierr; 2094 PetscInt i,j; 2095 2096 PetscFunctionBegin; 2097 if (type == NORM_FROBENIUS) { 2098 #if defined(PETSC_USE_REAL___FP16) 2099 PetscBLASInt one = 1,nz = a->nz; 2100 *nrm = BLASnrm2_(&nz,v,&one); 2101 #else 2102 for (i=0; i<a->nz; i++) { 2103 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2104 } 2105 *nrm = PetscSqrtReal(sum); 2106 #endif 2107 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 2108 } else if (type == NORM_1) { 2109 PetscReal *tmp; 2110 PetscInt *jj = a->j; 2111 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 2112 *nrm = 0.0; 2113 for (j=0; j<a->nz; j++) { 2114 tmp[*jj++] += PetscAbsScalar(*v); v++; 2115 } 2116 for (j=0; j<A->cmap->n; j++) { 2117 if (tmp[j] > *nrm) *nrm = tmp[j]; 2118 } 2119 ierr = PetscFree(tmp);CHKERRQ(ierr); 2120 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 2121 } else if (type == NORM_INFINITY) { 2122 *nrm = 0.0; 2123 for (j=0; j<A->rmap->n; j++) { 2124 v = a->a + a->i[j]; 2125 sum = 0.0; 2126 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2127 sum += PetscAbsScalar(*v); v++; 2128 } 2129 if (sum > *nrm) *nrm = sum; 2130 } 2131 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 2132 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2133 PetscFunctionReturn(0); 2134 } 2135 2136 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2137 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2138 { 2139 PetscErrorCode ierr; 2140 PetscInt i,j,anzj; 2141 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2142 PetscInt an=A->cmap->N,am=A->rmap->N; 2143 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2144 2145 PetscFunctionBegin; 2146 /* Allocate space for symbolic transpose info and work array */ 2147 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 2148 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 2149 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 2150 2151 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2152 /* Note: offset by 1 for fast conversion into csr format. */ 2153 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2154 /* Form ati for csr format of A^T. */ 2155 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2156 2157 /* Copy ati into atfill so we have locations of the next free space in atj */ 2158 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2159 2160 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2161 for (i=0;i<am;i++) { 2162 anzj = ai[i+1] - ai[i]; 2163 for (j=0;j<anzj;j++) { 2164 atj[atfill[*aj]] = i; 2165 atfill[*aj++] += 1; 2166 } 2167 } 2168 2169 /* Clean up temporary space and complete requests. */ 2170 ierr = PetscFree(atfill);CHKERRQ(ierr); 2171 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2172 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2173 2174 b = (Mat_SeqAIJ*)((*B)->data); 2175 b->free_a = PETSC_FALSE; 2176 b->free_ij = PETSC_TRUE; 2177 b->nonew = 0; 2178 PetscFunctionReturn(0); 2179 } 2180 2181 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2182 { 2183 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2184 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2185 MatScalar *va,*vb; 2186 PetscErrorCode ierr; 2187 PetscInt ma,na,mb,nb, i; 2188 2189 PetscFunctionBegin; 2190 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2191 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2192 if (ma!=nb || na!=mb) { 2193 *f = PETSC_FALSE; 2194 PetscFunctionReturn(0); 2195 } 2196 aii = aij->i; bii = bij->i; 2197 adx = aij->j; bdx = bij->j; 2198 va = aij->a; vb = bij->a; 2199 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2200 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2201 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2202 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2203 2204 *f = PETSC_TRUE; 2205 for (i=0; i<ma; i++) { 2206 while (aptr[i]<aii[i+1]) { 2207 PetscInt idc,idr; 2208 PetscScalar vc,vr; 2209 /* column/row index/value */ 2210 idc = adx[aptr[i]]; 2211 idr = bdx[bptr[idc]]; 2212 vc = va[aptr[i]]; 2213 vr = vb[bptr[idc]]; 2214 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2215 *f = PETSC_FALSE; 2216 goto done; 2217 } else { 2218 aptr[i]++; 2219 if (B || i!=idc) bptr[idc]++; 2220 } 2221 } 2222 } 2223 done: 2224 ierr = PetscFree(aptr);CHKERRQ(ierr); 2225 ierr = PetscFree(bptr);CHKERRQ(ierr); 2226 PetscFunctionReturn(0); 2227 } 2228 2229 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2230 { 2231 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2232 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2233 MatScalar *va,*vb; 2234 PetscErrorCode ierr; 2235 PetscInt ma,na,mb,nb, i; 2236 2237 PetscFunctionBegin; 2238 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2239 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2240 if (ma!=nb || na!=mb) { 2241 *f = PETSC_FALSE; 2242 PetscFunctionReturn(0); 2243 } 2244 aii = aij->i; bii = bij->i; 2245 adx = aij->j; bdx = bij->j; 2246 va = aij->a; vb = bij->a; 2247 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2248 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2249 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2250 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2251 2252 *f = PETSC_TRUE; 2253 for (i=0; i<ma; i++) { 2254 while (aptr[i]<aii[i+1]) { 2255 PetscInt idc,idr; 2256 PetscScalar vc,vr; 2257 /* column/row index/value */ 2258 idc = adx[aptr[i]]; 2259 idr = bdx[bptr[idc]]; 2260 vc = va[aptr[i]]; 2261 vr = vb[bptr[idc]]; 2262 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2263 *f = PETSC_FALSE; 2264 goto done; 2265 } else { 2266 aptr[i]++; 2267 if (B || i!=idc) bptr[idc]++; 2268 } 2269 } 2270 } 2271 done: 2272 ierr = PetscFree(aptr);CHKERRQ(ierr); 2273 ierr = PetscFree(bptr);CHKERRQ(ierr); 2274 PetscFunctionReturn(0); 2275 } 2276 2277 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2278 { 2279 PetscErrorCode ierr; 2280 2281 PetscFunctionBegin; 2282 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2283 PetscFunctionReturn(0); 2284 } 2285 2286 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2287 { 2288 PetscErrorCode ierr; 2289 2290 PetscFunctionBegin; 2291 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2292 PetscFunctionReturn(0); 2293 } 2294 2295 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2296 { 2297 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2298 const PetscScalar *l,*r; 2299 PetscScalar x; 2300 MatScalar *v; 2301 PetscErrorCode ierr; 2302 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2303 const PetscInt *jj; 2304 2305 PetscFunctionBegin; 2306 if (ll) { 2307 /* The local size is used so that VecMPI can be passed to this routine 2308 by MatDiagonalScale_MPIAIJ */ 2309 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2310 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2311 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 2312 v = a->a; 2313 for (i=0; i<m; i++) { 2314 x = l[i]; 2315 M = a->i[i+1] - a->i[i]; 2316 for (j=0; j<M; j++) (*v++) *= x; 2317 } 2318 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 2319 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2320 } 2321 if (rr) { 2322 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2323 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2324 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 2325 v = a->a; jj = a->j; 2326 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2327 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 2328 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2329 } 2330 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2331 PetscFunctionReturn(0); 2332 } 2333 2334 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2335 { 2336 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2337 PetscErrorCode ierr; 2338 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2339 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2340 const PetscInt *irow,*icol; 2341 PetscInt nrows,ncols; 2342 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2343 MatScalar *a_new,*mat_a; 2344 Mat C; 2345 PetscBool stride; 2346 2347 PetscFunctionBegin; 2348 2349 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2350 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2351 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2352 2353 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2354 if (stride) { 2355 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2356 } else { 2357 first = 0; 2358 step = 0; 2359 } 2360 if (stride && step == 1) { 2361 /* special case of contiguous rows */ 2362 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2363 /* loop over new rows determining lens and starting points */ 2364 for (i=0; i<nrows; i++) { 2365 kstart = ai[irow[i]]; 2366 kend = kstart + ailen[irow[i]]; 2367 starts[i] = kstart; 2368 for (k=kstart; k<kend; k++) { 2369 if (aj[k] >= first) { 2370 starts[i] = k; 2371 break; 2372 } 2373 } 2374 sum = 0; 2375 while (k < kend) { 2376 if (aj[k++] >= first+ncols) break; 2377 sum++; 2378 } 2379 lens[i] = sum; 2380 } 2381 /* create submatrix */ 2382 if (scall == MAT_REUSE_MATRIX) { 2383 PetscInt n_cols,n_rows; 2384 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2385 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2386 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2387 C = *B; 2388 } else { 2389 PetscInt rbs,cbs; 2390 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2391 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2392 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2393 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2394 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2395 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2396 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2397 } 2398 c = (Mat_SeqAIJ*)C->data; 2399 2400 /* loop over rows inserting into submatrix */ 2401 a_new = c->a; 2402 j_new = c->j; 2403 i_new = c->i; 2404 2405 for (i=0; i<nrows; i++) { 2406 ii = starts[i]; 2407 lensi = lens[i]; 2408 for (k=0; k<lensi; k++) { 2409 *j_new++ = aj[ii+k] - first; 2410 } 2411 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2412 a_new += lensi; 2413 i_new[i+1] = i_new[i] + lensi; 2414 c->ilen[i] = lensi; 2415 } 2416 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2417 } else { 2418 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2419 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2420 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2421 for (i=0; i<ncols; i++) { 2422 #if defined(PETSC_USE_DEBUG) 2423 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); 2424 #endif 2425 smap[icol[i]] = i+1; 2426 } 2427 2428 /* determine lens of each row */ 2429 for (i=0; i<nrows; i++) { 2430 kstart = ai[irow[i]]; 2431 kend = kstart + a->ilen[irow[i]]; 2432 lens[i] = 0; 2433 for (k=kstart; k<kend; k++) { 2434 if (smap[aj[k]]) { 2435 lens[i]++; 2436 } 2437 } 2438 } 2439 /* Create and fill new matrix */ 2440 if (scall == MAT_REUSE_MATRIX) { 2441 PetscBool equal; 2442 2443 c = (Mat_SeqAIJ*)((*B)->data); 2444 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2445 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2446 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2447 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2448 C = *B; 2449 } else { 2450 PetscInt rbs,cbs; 2451 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2452 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2453 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2454 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2455 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2456 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2457 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2458 } 2459 c = (Mat_SeqAIJ*)(C->data); 2460 for (i=0; i<nrows; i++) { 2461 row = irow[i]; 2462 kstart = ai[row]; 2463 kend = kstart + a->ilen[row]; 2464 mat_i = c->i[i]; 2465 mat_j = c->j + mat_i; 2466 mat_a = c->a + mat_i; 2467 mat_ilen = c->ilen + i; 2468 for (k=kstart; k<kend; k++) { 2469 if ((tcol=smap[a->j[k]])) { 2470 *mat_j++ = tcol - 1; 2471 *mat_a++ = a->a[k]; 2472 (*mat_ilen)++; 2473 2474 } 2475 } 2476 } 2477 /* Free work space */ 2478 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2479 ierr = PetscFree(smap);CHKERRQ(ierr); 2480 ierr = PetscFree(lens);CHKERRQ(ierr); 2481 /* sort */ 2482 for (i = 0; i < nrows; i++) { 2483 PetscInt ilen; 2484 2485 mat_i = c->i[i]; 2486 mat_j = c->j + mat_i; 2487 mat_a = c->a + mat_i; 2488 ilen = c->ilen[i]; 2489 ierr = PetscSortIntWithScalarArray(ilen,mat_j,mat_a);CHKERRQ(ierr); 2490 } 2491 } 2492 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2493 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2494 2495 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2496 *B = C; 2497 PetscFunctionReturn(0); 2498 } 2499 2500 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2501 { 2502 PetscErrorCode ierr; 2503 Mat B; 2504 2505 PetscFunctionBegin; 2506 if (scall == MAT_INITIAL_MATRIX) { 2507 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2508 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2509 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2510 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2511 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2512 *subMat = B; 2513 } else { 2514 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2515 } 2516 PetscFunctionReturn(0); 2517 } 2518 2519 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2520 { 2521 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2522 PetscErrorCode ierr; 2523 Mat outA; 2524 PetscBool row_identity,col_identity; 2525 2526 PetscFunctionBegin; 2527 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2528 2529 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2530 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2531 2532 outA = inA; 2533 outA->factortype = MAT_FACTOR_LU; 2534 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 2535 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 2536 2537 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2538 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2539 2540 a->row = row; 2541 2542 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2543 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2544 2545 a->col = col; 2546 2547 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2548 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2549 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2550 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2551 2552 if (!a->solve_work) { /* this matrix may have been factored before */ 2553 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2554 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2555 } 2556 2557 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2558 if (row_identity && col_identity) { 2559 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2560 } else { 2561 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2562 } 2563 PetscFunctionReturn(0); 2564 } 2565 2566 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2567 { 2568 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2569 PetscScalar oalpha = alpha; 2570 PetscErrorCode ierr; 2571 PetscBLASInt one = 1,bnz; 2572 2573 PetscFunctionBegin; 2574 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2575 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2576 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2577 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2578 PetscFunctionReturn(0); 2579 } 2580 2581 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2582 { 2583 PetscErrorCode ierr; 2584 PetscInt i; 2585 2586 PetscFunctionBegin; 2587 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2588 ierr = PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr);CHKERRQ(ierr); 2589 2590 for (i=0; i<submatj->nrqr; ++i) { 2591 ierr = PetscFree(submatj->sbuf2[i]);CHKERRQ(ierr); 2592 } 2593 ierr = PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1);CHKERRQ(ierr); 2594 2595 if (submatj->rbuf1) { 2596 ierr = PetscFree(submatj->rbuf1[0]);CHKERRQ(ierr); 2597 ierr = PetscFree(submatj->rbuf1);CHKERRQ(ierr); 2598 } 2599 2600 for (i=0; i<submatj->nrqs; ++i) { 2601 ierr = PetscFree(submatj->rbuf3[i]);CHKERRQ(ierr); 2602 } 2603 ierr = PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3);CHKERRQ(ierr); 2604 ierr = PetscFree(submatj->pa);CHKERRQ(ierr); 2605 } 2606 2607 #if defined(PETSC_USE_CTABLE) 2608 ierr = PetscTableDestroy((PetscTable*)&submatj->rmap);CHKERRQ(ierr); 2609 if (submatj->cmap_loc) {ierr = PetscFree(submatj->cmap_loc);CHKERRQ(ierr);} 2610 ierr = PetscFree(submatj->rmap_loc);CHKERRQ(ierr); 2611 #else 2612 ierr = PetscFree(submatj->rmap);CHKERRQ(ierr); 2613 #endif 2614 2615 if (!submatj->allcolumns) { 2616 #if defined(PETSC_USE_CTABLE) 2617 ierr = PetscTableDestroy((PetscTable*)&submatj->cmap);CHKERRQ(ierr); 2618 #else 2619 ierr = PetscFree(submatj->cmap);CHKERRQ(ierr); 2620 #endif 2621 } 2622 ierr = PetscFree(submatj->row2proc);CHKERRQ(ierr); 2623 2624 ierr = PetscFree(submatj);CHKERRQ(ierr); 2625 PetscFunctionReturn(0); 2626 } 2627 2628 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2629 { 2630 PetscErrorCode ierr; 2631 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2632 Mat_SubSppt *submatj = c->submatis1; 2633 2634 PetscFunctionBegin; 2635 ierr = (*submatj->destroy)(C);CHKERRQ(ierr); 2636 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2637 PetscFunctionReturn(0); 2638 } 2639 2640 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2641 { 2642 PetscErrorCode ierr; 2643 PetscInt i; 2644 Mat C; 2645 Mat_SeqAIJ *c; 2646 Mat_SubSppt *submatj; 2647 2648 PetscFunctionBegin; 2649 for (i=0; i<n; i++) { 2650 C = (*mat)[i]; 2651 c = (Mat_SeqAIJ*)C->data; 2652 submatj = c->submatis1; 2653 if (submatj) { 2654 if (--((PetscObject)C)->refct <= 0) { 2655 ierr = (*submatj->destroy)(C);CHKERRQ(ierr); 2656 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2657 ierr = PetscFree(C->defaultvectype);CHKERRQ(ierr); 2658 ierr = PetscLayoutDestroy(&C->rmap);CHKERRQ(ierr); 2659 ierr = PetscLayoutDestroy(&C->cmap);CHKERRQ(ierr); 2660 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 2661 } 2662 } else { 2663 ierr = MatDestroy(&C);CHKERRQ(ierr); 2664 } 2665 } 2666 2667 /* Destroy Dummy submatrices created for reuse */ 2668 ierr = MatDestroySubMatrices_Dummy(n,mat);CHKERRQ(ierr); 2669 2670 ierr = PetscFree(*mat);CHKERRQ(ierr); 2671 PetscFunctionReturn(0); 2672 } 2673 2674 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2675 { 2676 PetscErrorCode ierr; 2677 PetscInt i; 2678 2679 PetscFunctionBegin; 2680 if (scall == MAT_INITIAL_MATRIX) { 2681 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2682 } 2683 2684 for (i=0; i<n; i++) { 2685 ierr = MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2686 } 2687 PetscFunctionReturn(0); 2688 } 2689 2690 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2691 { 2692 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2693 PetscErrorCode ierr; 2694 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2695 const PetscInt *idx; 2696 PetscInt start,end,*ai,*aj; 2697 PetscBT table; 2698 2699 PetscFunctionBegin; 2700 m = A->rmap->n; 2701 ai = a->i; 2702 aj = a->j; 2703 2704 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2705 2706 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2707 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2708 2709 for (i=0; i<is_max; i++) { 2710 /* Initialize the two local arrays */ 2711 isz = 0; 2712 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2713 2714 /* Extract the indices, assume there can be duplicate entries */ 2715 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2716 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2717 2718 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2719 for (j=0; j<n; ++j) { 2720 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2721 } 2722 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2723 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2724 2725 k = 0; 2726 for (j=0; j<ov; j++) { /* for each overlap */ 2727 n = isz; 2728 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2729 row = nidx[k]; 2730 start = ai[row]; 2731 end = ai[row+1]; 2732 for (l = start; l<end; l++) { 2733 val = aj[l]; 2734 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2735 } 2736 } 2737 } 2738 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2739 } 2740 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2741 ierr = PetscFree(nidx);CHKERRQ(ierr); 2742 PetscFunctionReturn(0); 2743 } 2744 2745 /* -------------------------------------------------------------- */ 2746 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2747 { 2748 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2749 PetscErrorCode ierr; 2750 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2751 const PetscInt *row,*col; 2752 PetscInt *cnew,j,*lens; 2753 IS icolp,irowp; 2754 PetscInt *cwork = NULL; 2755 PetscScalar *vwork = NULL; 2756 2757 PetscFunctionBegin; 2758 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2759 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2760 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2761 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2762 2763 /* determine lengths of permuted rows */ 2764 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2765 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2766 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2767 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2768 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2769 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2770 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2771 ierr = PetscFree(lens);CHKERRQ(ierr); 2772 2773 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2774 for (i=0; i<m; i++) { 2775 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2776 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2777 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2778 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2779 } 2780 ierr = PetscFree(cnew);CHKERRQ(ierr); 2781 2782 (*B)->assembled = PETSC_FALSE; 2783 2784 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2785 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2786 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2787 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2788 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2789 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2790 PetscFunctionReturn(0); 2791 } 2792 2793 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2794 { 2795 PetscErrorCode ierr; 2796 2797 PetscFunctionBegin; 2798 /* If the two matrices have the same copy implementation, use fast copy. */ 2799 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2800 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2801 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2802 2803 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"); 2804 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2805 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2806 } else { 2807 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2808 } 2809 PetscFunctionReturn(0); 2810 } 2811 2812 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2813 { 2814 PetscErrorCode ierr; 2815 2816 PetscFunctionBegin; 2817 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2818 PetscFunctionReturn(0); 2819 } 2820 2821 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2822 { 2823 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2824 2825 PetscFunctionBegin; 2826 *array = a->a; 2827 PetscFunctionReturn(0); 2828 } 2829 2830 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2831 { 2832 PetscFunctionBegin; 2833 PetscFunctionReturn(0); 2834 } 2835 2836 /* 2837 Computes the number of nonzeros per row needed for preallocation when X and Y 2838 have different nonzero structure. 2839 */ 2840 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2841 { 2842 PetscInt i,j,k,nzx,nzy; 2843 2844 PetscFunctionBegin; 2845 /* Set the number of nonzeros in the new matrix */ 2846 for (i=0; i<m; i++) { 2847 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2848 nzx = xi[i+1] - xi[i]; 2849 nzy = yi[i+1] - yi[i]; 2850 nnz[i] = 0; 2851 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2852 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2853 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2854 nnz[i]++; 2855 } 2856 for (; k<nzy; k++) nnz[i]++; 2857 } 2858 PetscFunctionReturn(0); 2859 } 2860 2861 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2862 { 2863 PetscInt m = Y->rmap->N; 2864 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2865 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2866 PetscErrorCode ierr; 2867 2868 PetscFunctionBegin; 2869 /* Set the number of nonzeros in the new matrix */ 2870 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2871 PetscFunctionReturn(0); 2872 } 2873 2874 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2875 { 2876 PetscErrorCode ierr; 2877 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2878 PetscBLASInt one=1,bnz; 2879 2880 PetscFunctionBegin; 2881 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2882 if (str == SAME_NONZERO_PATTERN) { 2883 PetscScalar alpha = a; 2884 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2885 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2886 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2887 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2888 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2889 } else { 2890 Mat B; 2891 PetscInt *nnz; 2892 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2893 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2894 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2895 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2896 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2897 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2898 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2899 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2900 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2901 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2902 ierr = PetscFree(nnz);CHKERRQ(ierr); 2903 } 2904 PetscFunctionReturn(0); 2905 } 2906 2907 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2908 { 2909 #if defined(PETSC_USE_COMPLEX) 2910 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2911 PetscInt i,nz; 2912 PetscScalar *a; 2913 2914 PetscFunctionBegin; 2915 nz = aij->nz; 2916 a = aij->a; 2917 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2918 #else 2919 PetscFunctionBegin; 2920 #endif 2921 PetscFunctionReturn(0); 2922 } 2923 2924 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2925 { 2926 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2927 PetscErrorCode ierr; 2928 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2929 PetscReal atmp; 2930 PetscScalar *x; 2931 MatScalar *aa; 2932 2933 PetscFunctionBegin; 2934 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2935 aa = a->a; 2936 ai = a->i; 2937 aj = a->j; 2938 2939 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2940 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2941 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2942 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2943 for (i=0; i<m; i++) { 2944 ncols = ai[1] - ai[0]; ai++; 2945 x[i] = 0.0; 2946 for (j=0; j<ncols; j++) { 2947 atmp = PetscAbsScalar(*aa); 2948 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2949 aa++; aj++; 2950 } 2951 } 2952 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2953 PetscFunctionReturn(0); 2954 } 2955 2956 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2957 { 2958 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2959 PetscErrorCode ierr; 2960 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2961 PetscScalar *x; 2962 MatScalar *aa; 2963 2964 PetscFunctionBegin; 2965 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2966 aa = a->a; 2967 ai = a->i; 2968 aj = a->j; 2969 2970 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2971 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2972 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2973 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2974 for (i=0; i<m; i++) { 2975 ncols = ai[1] - ai[0]; ai++; 2976 if (ncols == A->cmap->n) { /* row is dense */ 2977 x[i] = *aa; if (idx) idx[i] = 0; 2978 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2979 x[i] = 0.0; 2980 if (idx) { 2981 idx[i] = 0; /* in case ncols is zero */ 2982 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2983 if (aj[j] > j) { 2984 idx[i] = j; 2985 break; 2986 } 2987 } 2988 } 2989 } 2990 for (j=0; j<ncols; j++) { 2991 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2992 aa++; aj++; 2993 } 2994 } 2995 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2996 PetscFunctionReturn(0); 2997 } 2998 2999 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3000 { 3001 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3002 PetscErrorCode ierr; 3003 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3004 PetscReal atmp; 3005 PetscScalar *x; 3006 MatScalar *aa; 3007 3008 PetscFunctionBegin; 3009 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3010 aa = a->a; 3011 ai = a->i; 3012 aj = a->j; 3013 3014 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3015 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3016 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3017 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); 3018 for (i=0; i<m; i++) { 3019 ncols = ai[1] - ai[0]; ai++; 3020 if (ncols) { 3021 /* Get first nonzero */ 3022 for (j = 0; j < ncols; j++) { 3023 atmp = PetscAbsScalar(aa[j]); 3024 if (atmp > 1.0e-12) { 3025 x[i] = atmp; 3026 if (idx) idx[i] = aj[j]; 3027 break; 3028 } 3029 } 3030 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 3031 } else { 3032 x[i] = 0.0; if (idx) idx[i] = 0; 3033 } 3034 for (j = 0; j < ncols; j++) { 3035 atmp = PetscAbsScalar(*aa); 3036 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3037 aa++; aj++; 3038 } 3039 } 3040 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3041 PetscFunctionReturn(0); 3042 } 3043 3044 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3045 { 3046 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3047 PetscErrorCode ierr; 3048 PetscInt i,j,m = A->rmap->n,ncols,n; 3049 const PetscInt *ai,*aj; 3050 PetscScalar *x; 3051 const MatScalar *aa; 3052 3053 PetscFunctionBegin; 3054 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3055 aa = a->a; 3056 ai = a->i; 3057 aj = a->j; 3058 3059 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3060 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3061 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3062 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3063 for (i=0; i<m; i++) { 3064 ncols = ai[1] - ai[0]; ai++; 3065 if (ncols == A->cmap->n) { /* row is dense */ 3066 x[i] = *aa; if (idx) idx[i] = 0; 3067 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3068 x[i] = 0.0; 3069 if (idx) { /* find first implicit 0.0 in the row */ 3070 idx[i] = 0; /* in case ncols is zero */ 3071 for (j=0; j<ncols; j++) { 3072 if (aj[j] > j) { 3073 idx[i] = j; 3074 break; 3075 } 3076 } 3077 } 3078 } 3079 for (j=0; j<ncols; j++) { 3080 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3081 aa++; aj++; 3082 } 3083 } 3084 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3085 PetscFunctionReturn(0); 3086 } 3087 3088 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3089 { 3090 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3091 PetscErrorCode ierr; 3092 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3093 MatScalar *diag,work[25],*v_work; 3094 const PetscReal shift = 0.0; 3095 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3096 3097 PetscFunctionBegin; 3098 allowzeropivot = PetscNot(A->erroriffailure); 3099 if (a->ibdiagvalid) { 3100 if (values) *values = a->ibdiag; 3101 PetscFunctionReturn(0); 3102 } 3103 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 3104 if (!a->ibdiag) { 3105 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 3106 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 3107 } 3108 diag = a->ibdiag; 3109 if (values) *values = a->ibdiag; 3110 /* factor and invert each block */ 3111 switch (bs) { 3112 case 1: 3113 for (i=0; i<mbs; i++) { 3114 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3115 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3116 if (allowzeropivot) { 3117 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3118 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3119 A->factorerror_zeropivot_row = i; 3120 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 3121 } 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); 3122 } 3123 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3124 } 3125 break; 3126 case 2: 3127 for (i=0; i<mbs; i++) { 3128 ij[0] = 2*i; ij[1] = 2*i + 1; 3129 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3130 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3131 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3132 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3133 diag += 4; 3134 } 3135 break; 3136 case 3: 3137 for (i=0; i<mbs; i++) { 3138 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3139 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3140 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3141 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3142 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3143 diag += 9; 3144 } 3145 break; 3146 case 4: 3147 for (i=0; i<mbs; i++) { 3148 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3149 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3150 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3151 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3152 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3153 diag += 16; 3154 } 3155 break; 3156 case 5: 3157 for (i=0; i<mbs; i++) { 3158 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3159 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3160 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3161 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3162 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3163 diag += 25; 3164 } 3165 break; 3166 case 6: 3167 for (i=0; i<mbs; i++) { 3168 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; 3169 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3170 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3171 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3172 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3173 diag += 36; 3174 } 3175 break; 3176 case 7: 3177 for (i=0; i<mbs; i++) { 3178 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; 3179 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3180 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3181 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3182 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3183 diag += 49; 3184 } 3185 break; 3186 default: 3187 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3188 for (i=0; i<mbs; i++) { 3189 for (j=0; j<bs; j++) { 3190 IJ[j] = bs*i + j; 3191 } 3192 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3193 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3194 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3195 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3196 diag += bs2; 3197 } 3198 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3199 } 3200 a->ibdiagvalid = PETSC_TRUE; 3201 PetscFunctionReturn(0); 3202 } 3203 3204 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3205 { 3206 PetscErrorCode ierr; 3207 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3208 PetscScalar a; 3209 PetscInt m,n,i,j,col; 3210 3211 PetscFunctionBegin; 3212 if (!x->assembled) { 3213 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3214 for (i=0; i<m; i++) { 3215 for (j=0; j<aij->imax[i]; j++) { 3216 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3217 col = (PetscInt)(n*PetscRealPart(a)); 3218 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3219 } 3220 } 3221 } else { 3222 for (i=0; i<aij->nz; i++) {ierr = PetscRandomGetValue(rctx,aij->a+i);CHKERRQ(ierr);} 3223 } 3224 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3225 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3226 PetscFunctionReturn(0); 3227 } 3228 3229 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3230 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx) 3231 { 3232 PetscErrorCode ierr; 3233 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3234 PetscScalar a; 3235 PetscInt m,n,i,j,col,nskip; 3236 3237 PetscFunctionBegin; 3238 nskip = high - low; 3239 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3240 n -= nskip; /* shrink number of columns where nonzeros can be set */ 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 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3246 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3247 } 3248 } 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 /* -------------------------------------------------------------------*/ 3256 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3257 MatGetRow_SeqAIJ, 3258 MatRestoreRow_SeqAIJ, 3259 MatMult_SeqAIJ, 3260 /* 4*/ MatMultAdd_SeqAIJ, 3261 MatMultTranspose_SeqAIJ, 3262 MatMultTransposeAdd_SeqAIJ, 3263 0, 3264 0, 3265 0, 3266 /* 10*/ 0, 3267 MatLUFactor_SeqAIJ, 3268 0, 3269 MatSOR_SeqAIJ, 3270 MatTranspose_SeqAIJ, 3271 /*1 5*/ MatGetInfo_SeqAIJ, 3272 MatEqual_SeqAIJ, 3273 MatGetDiagonal_SeqAIJ, 3274 MatDiagonalScale_SeqAIJ, 3275 MatNorm_SeqAIJ, 3276 /* 20*/ 0, 3277 MatAssemblyEnd_SeqAIJ, 3278 MatSetOption_SeqAIJ, 3279 MatZeroEntries_SeqAIJ, 3280 /* 24*/ MatZeroRows_SeqAIJ, 3281 0, 3282 0, 3283 0, 3284 0, 3285 /* 29*/ MatSetUp_SeqAIJ, 3286 0, 3287 0, 3288 0, 3289 0, 3290 /* 34*/ MatDuplicate_SeqAIJ, 3291 0, 3292 0, 3293 MatILUFactor_SeqAIJ, 3294 0, 3295 /* 39*/ MatAXPY_SeqAIJ, 3296 MatCreateSubMatrices_SeqAIJ, 3297 MatIncreaseOverlap_SeqAIJ, 3298 MatGetValues_SeqAIJ, 3299 MatCopy_SeqAIJ, 3300 /* 44*/ MatGetRowMax_SeqAIJ, 3301 MatScale_SeqAIJ, 3302 MatShift_SeqAIJ, 3303 MatDiagonalSet_SeqAIJ, 3304 MatZeroRowsColumns_SeqAIJ, 3305 /* 49*/ MatSetRandom_SeqAIJ, 3306 MatGetRowIJ_SeqAIJ, 3307 MatRestoreRowIJ_SeqAIJ, 3308 MatGetColumnIJ_SeqAIJ, 3309 MatRestoreColumnIJ_SeqAIJ, 3310 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3311 0, 3312 0, 3313 MatPermute_SeqAIJ, 3314 0, 3315 /* 59*/ 0, 3316 MatDestroy_SeqAIJ, 3317 MatView_SeqAIJ, 3318 0, 3319 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3320 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3321 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3322 0, 3323 0, 3324 0, 3325 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3326 MatGetRowMinAbs_SeqAIJ, 3327 0, 3328 0, 3329 0, 3330 /* 74*/ 0, 3331 MatFDColoringApply_AIJ, 3332 0, 3333 0, 3334 0, 3335 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3336 0, 3337 0, 3338 0, 3339 MatLoad_SeqAIJ, 3340 /* 84*/ MatIsSymmetric_SeqAIJ, 3341 MatIsHermitian_SeqAIJ, 3342 0, 3343 0, 3344 0, 3345 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3346 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3347 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3348 MatPtAP_SeqAIJ_SeqAIJ, 3349 MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy, 3350 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3351 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3352 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3353 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3354 0, 3355 /* 99*/ 0, 3356 0, 3357 0, 3358 MatConjugate_SeqAIJ, 3359 0, 3360 /*104*/ MatSetValuesRow_SeqAIJ, 3361 MatRealPart_SeqAIJ, 3362 MatImaginaryPart_SeqAIJ, 3363 0, 3364 0, 3365 /*109*/ MatMatSolve_SeqAIJ, 3366 0, 3367 MatGetRowMin_SeqAIJ, 3368 0, 3369 MatMissingDiagonal_SeqAIJ, 3370 /*114*/ 0, 3371 0, 3372 0, 3373 0, 3374 0, 3375 /*119*/ 0, 3376 0, 3377 0, 3378 0, 3379 MatGetMultiProcBlock_SeqAIJ, 3380 /*124*/ MatFindNonzeroRows_SeqAIJ, 3381 MatGetColumnNorms_SeqAIJ, 3382 MatInvertBlockDiagonal_SeqAIJ, 3383 MatInvertVariableBlockDiagonal_SeqAIJ, 3384 0, 3385 /*129*/ 0, 3386 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3387 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3388 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3389 MatTransposeColoringCreate_SeqAIJ, 3390 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3391 MatTransColoringApplyDenToSp_SeqAIJ, 3392 MatRARt_SeqAIJ_SeqAIJ, 3393 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3394 MatRARtNumeric_SeqAIJ_SeqAIJ, 3395 /*139*/0, 3396 0, 3397 0, 3398 MatFDColoringSetUp_SeqXAIJ, 3399 MatFindOffBlockDiagonalEntries_SeqAIJ, 3400 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3401 MatDestroySubMatrices_SeqAIJ 3402 }; 3403 3404 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3405 { 3406 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3407 PetscInt i,nz,n; 3408 3409 PetscFunctionBegin; 3410 nz = aij->maxnz; 3411 n = mat->rmap->n; 3412 for (i=0; i<nz; i++) { 3413 aij->j[i] = indices[i]; 3414 } 3415 aij->nz = nz; 3416 for (i=0; i<n; i++) { 3417 aij->ilen[i] = aij->imax[i]; 3418 } 3419 PetscFunctionReturn(0); 3420 } 3421 3422 /* 3423 * When a sparse matrix has many zero columns, we should compact them out to save the space 3424 * This happens in MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3425 * */ 3426 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3427 { 3428 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3429 PetscTable gid1_lid1; 3430 PetscTablePosition tpos; 3431 PetscInt gid,lid,i,j,ncols,ec; 3432 PetscInt *garray; 3433 PetscErrorCode ierr; 3434 3435 PetscFunctionBegin; 3436 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3437 PetscValidPointer(mapping,2); 3438 /* use a table */ 3439 ierr = PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1);CHKERRQ(ierr); 3440 ec = 0; 3441 for (i=0; i<mat->rmap->n; i++) { 3442 ncols = aij->i[i+1] - aij->i[i]; 3443 for (j=0; j<ncols; j++) { 3444 PetscInt data,gid1 = aij->j[aij->i[i] + j] + 1; 3445 ierr = PetscTableFind(gid1_lid1,gid1,&data);CHKERRQ(ierr); 3446 if (!data) { 3447 /* one based table */ 3448 ierr = PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);CHKERRQ(ierr); 3449 } 3450 } 3451 } 3452 /* form array of columns we need */ 3453 ierr = PetscMalloc1(ec+1,&garray);CHKERRQ(ierr); 3454 ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr); 3455 while (tpos) { 3456 ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);CHKERRQ(ierr); 3457 gid--; 3458 lid--; 3459 garray[lid] = gid; 3460 } 3461 ierr = PetscSortInt(ec,garray);CHKERRQ(ierr); /* sort, and rebuild */ 3462 ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr); 3463 for (i=0; i<ec; i++) { 3464 ierr = PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);CHKERRQ(ierr); 3465 } 3466 /* compact out the extra columns in B */ 3467 for (i=0; i<mat->rmap->n; i++) { 3468 ncols = aij->i[i+1] - aij->i[i]; 3469 for (j=0; j<ncols; j++) { 3470 PetscInt gid1 = aij->j[aij->i[i] + j] + 1; 3471 ierr = PetscTableFind(gid1_lid1,gid1,&lid);CHKERRQ(ierr); 3472 lid--; 3473 aij->j[aij->i[i] + j] = lid; 3474 } 3475 } 3476 mat->cmap->n = mat->cmap->N = ec; 3477 mat->cmap->bs = 1; 3478 3479 ierr = PetscTableDestroy(&gid1_lid1);CHKERRQ(ierr); 3480 ierr = PetscLayoutSetUp((mat->cmap));CHKERRQ(ierr); 3481 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping);CHKERRQ(ierr); 3482 ierr = ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH);CHKERRQ(ierr); 3483 PetscFunctionReturn(0); 3484 } 3485 3486 /*@ 3487 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3488 in the matrix. 3489 3490 Input Parameters: 3491 + mat - the SeqAIJ matrix 3492 - indices - the column indices 3493 3494 Level: advanced 3495 3496 Notes: 3497 This can be called if you have precomputed the nonzero structure of the 3498 matrix and want to provide it to the matrix object to improve the performance 3499 of the MatSetValues() operation. 3500 3501 You MUST have set the correct numbers of nonzeros per row in the call to 3502 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3503 3504 MUST be called before any calls to MatSetValues(); 3505 3506 The indices should start with zero, not one. 3507 3508 @*/ 3509 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3510 { 3511 PetscErrorCode ierr; 3512 3513 PetscFunctionBegin; 3514 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3515 PetscValidPointer(indices,2); 3516 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3517 PetscFunctionReturn(0); 3518 } 3519 3520 /* ----------------------------------------------------------------------------------------*/ 3521 3522 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3523 { 3524 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3525 PetscErrorCode ierr; 3526 size_t nz = aij->i[mat->rmap->n]; 3527 3528 PetscFunctionBegin; 3529 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3530 3531 /* allocate space for values if not already there */ 3532 if (!aij->saved_values) { 3533 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3534 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3535 } 3536 3537 /* copy values over */ 3538 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3539 PetscFunctionReturn(0); 3540 } 3541 3542 /*@ 3543 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3544 example, reuse of the linear part of a Jacobian, while recomputing the 3545 nonlinear portion. 3546 3547 Collect on Mat 3548 3549 Input Parameters: 3550 . mat - the matrix (currently only AIJ matrices support this option) 3551 3552 Level: advanced 3553 3554 Common Usage, with SNESSolve(): 3555 $ Create Jacobian matrix 3556 $ Set linear terms into matrix 3557 $ Apply boundary conditions to matrix, at this time matrix must have 3558 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3559 $ boundary conditions again will not change the nonzero structure 3560 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3561 $ ierr = MatStoreValues(mat); 3562 $ Call SNESSetJacobian() with matrix 3563 $ In your Jacobian routine 3564 $ ierr = MatRetrieveValues(mat); 3565 $ Set nonlinear terms in matrix 3566 3567 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3568 $ // build linear portion of Jacobian 3569 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3570 $ ierr = MatStoreValues(mat); 3571 $ loop over nonlinear iterations 3572 $ ierr = MatRetrieveValues(mat); 3573 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3574 $ // call MatAssemblyBegin/End() on matrix 3575 $ Solve linear system with Jacobian 3576 $ endloop 3577 3578 Notes: 3579 Matrix must already be assemblied before calling this routine 3580 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3581 calling this routine. 3582 3583 When this is called multiple times it overwrites the previous set of stored values 3584 and does not allocated additional space. 3585 3586 .seealso: MatRetrieveValues() 3587 3588 @*/ 3589 PetscErrorCode MatStoreValues(Mat mat) 3590 { 3591 PetscErrorCode ierr; 3592 3593 PetscFunctionBegin; 3594 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3595 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3596 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3597 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3598 PetscFunctionReturn(0); 3599 } 3600 3601 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3602 { 3603 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3604 PetscErrorCode ierr; 3605 PetscInt nz = aij->i[mat->rmap->n]; 3606 3607 PetscFunctionBegin; 3608 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3609 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3610 /* copy values over */ 3611 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3612 PetscFunctionReturn(0); 3613 } 3614 3615 /*@ 3616 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3617 example, reuse of the linear part of a Jacobian, while recomputing the 3618 nonlinear portion. 3619 3620 Collect on Mat 3621 3622 Input Parameters: 3623 . mat - the matrix (currently only AIJ matrices support this option) 3624 3625 Level: advanced 3626 3627 .seealso: MatStoreValues() 3628 3629 @*/ 3630 PetscErrorCode MatRetrieveValues(Mat mat) 3631 { 3632 PetscErrorCode ierr; 3633 3634 PetscFunctionBegin; 3635 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3636 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3637 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3638 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3639 PetscFunctionReturn(0); 3640 } 3641 3642 3643 /* --------------------------------------------------------------------------------*/ 3644 /*@C 3645 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3646 (the default parallel PETSc format). For good matrix assembly performance 3647 the user should preallocate the matrix storage by setting the parameter nz 3648 (or the array nnz). By setting these parameters accurately, performance 3649 during matrix assembly can be increased by more than a factor of 50. 3650 3651 Collective on MPI_Comm 3652 3653 Input Parameters: 3654 + comm - MPI communicator, set to PETSC_COMM_SELF 3655 . m - number of rows 3656 . n - number of columns 3657 . nz - number of nonzeros per row (same for all rows) 3658 - nnz - array containing the number of nonzeros in the various rows 3659 (possibly different for each row) or NULL 3660 3661 Output Parameter: 3662 . A - the matrix 3663 3664 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3665 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3666 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3667 3668 Notes: 3669 If nnz is given then nz is ignored 3670 3671 The AIJ format (also called the Yale sparse matrix format or 3672 compressed row storage), is fully compatible with standard Fortran 77 3673 storage. That is, the stored row and column indices can begin at 3674 either one (as in Fortran) or zero. See the users' manual for details. 3675 3676 Specify the preallocated storage with either nz or nnz (not both). 3677 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3678 allocation. For large problems you MUST preallocate memory or you 3679 will get TERRIBLE performance, see the users' manual chapter on matrices. 3680 3681 By default, this format uses inodes (identical nodes) when possible, to 3682 improve numerical efficiency of matrix-vector products and solves. We 3683 search for consecutive rows with the same nonzero structure, thereby 3684 reusing matrix information to achieve increased efficiency. 3685 3686 Options Database Keys: 3687 + -mat_no_inode - Do not use inodes 3688 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3689 3690 Level: intermediate 3691 3692 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3693 3694 @*/ 3695 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3696 { 3697 PetscErrorCode ierr; 3698 3699 PetscFunctionBegin; 3700 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3701 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3702 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3703 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3704 PetscFunctionReturn(0); 3705 } 3706 3707 /*@C 3708 MatSeqAIJSetPreallocation - For good matrix assembly performance 3709 the user should preallocate the matrix storage by setting the parameter nz 3710 (or the array nnz). By setting these parameters accurately, performance 3711 during matrix assembly can be increased by more than a factor of 50. 3712 3713 Collective on MPI_Comm 3714 3715 Input Parameters: 3716 + B - The matrix 3717 . nz - number of nonzeros per row (same for all rows) 3718 - nnz - array containing the number of nonzeros in the various rows 3719 (possibly different for each row) or NULL 3720 3721 Notes: 3722 If nnz is given then nz is ignored 3723 3724 The AIJ format (also called the Yale sparse matrix format or 3725 compressed row storage), is fully compatible with standard Fortran 77 3726 storage. That is, the stored row and column indices can begin at 3727 either one (as in Fortran) or zero. See the users' manual for details. 3728 3729 Specify the preallocated storage with either nz or nnz (not both). 3730 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3731 allocation. For large problems you MUST preallocate memory or you 3732 will get TERRIBLE performance, see the users' manual chapter on matrices. 3733 3734 You can call MatGetInfo() to get information on how effective the preallocation was; 3735 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3736 You can also run with the option -info and look for messages with the string 3737 malloc in them to see if additional memory allocation was needed. 3738 3739 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3740 entries or columns indices 3741 3742 By default, this format uses inodes (identical nodes) when possible, to 3743 improve numerical efficiency of matrix-vector products and solves. We 3744 search for consecutive rows with the same nonzero structure, thereby 3745 reusing matrix information to achieve increased efficiency. 3746 3747 Options Database Keys: 3748 + -mat_no_inode - Do not use inodes 3749 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3750 3751 Level: intermediate 3752 3753 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3754 3755 @*/ 3756 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3757 { 3758 PetscErrorCode ierr; 3759 3760 PetscFunctionBegin; 3761 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3762 PetscValidType(B,1); 3763 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3764 PetscFunctionReturn(0); 3765 } 3766 3767 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3768 { 3769 Mat_SeqAIJ *b; 3770 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3771 PetscErrorCode ierr; 3772 PetscInt i; 3773 3774 PetscFunctionBegin; 3775 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3776 if (nz == MAT_SKIP_ALLOCATION) { 3777 skipallocation = PETSC_TRUE; 3778 nz = 0; 3779 } 3780 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3781 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3782 3783 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3784 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3785 if (nnz) { 3786 for (i=0; i<B->rmap->n; i++) { 3787 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]); 3788 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); 3789 } 3790 } 3791 3792 B->preallocated = PETSC_TRUE; 3793 3794 b = (Mat_SeqAIJ*)B->data; 3795 3796 if (!skipallocation) { 3797 if (!b->imax) { 3798 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3799 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3800 } 3801 if (!b->ipre) { 3802 ierr = PetscMalloc1(B->rmap->n,&b->ipre);CHKERRQ(ierr); 3803 ierr = PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3804 } 3805 if (!nnz) { 3806 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3807 else if (nz < 0) nz = 1; 3808 nz = PetscMin(nz,B->cmap->n); 3809 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3810 nz = nz*B->rmap->n; 3811 } else { 3812 nz = 0; 3813 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3814 } 3815 /* b->ilen will count nonzeros in each row so far. */ 3816 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3817 3818 /* allocate the matrix space */ 3819 /* FIXME: should B's old memory be unlogged? */ 3820 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3821 if (B->structure_only) { 3822 ierr = PetscMalloc1(nz,&b->j);CHKERRQ(ierr); 3823 ierr = PetscMalloc1(B->rmap->n+1,&b->i);CHKERRQ(ierr); 3824 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));CHKERRQ(ierr); 3825 } else { 3826 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3827 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3828 } 3829 b->i[0] = 0; 3830 for (i=1; i<B->rmap->n+1; i++) { 3831 b->i[i] = b->i[i-1] + b->imax[i-1]; 3832 } 3833 if (B->structure_only) { 3834 b->singlemalloc = PETSC_FALSE; 3835 b->free_a = PETSC_FALSE; 3836 } else { 3837 b->singlemalloc = PETSC_TRUE; 3838 b->free_a = PETSC_TRUE; 3839 } 3840 b->free_ij = PETSC_TRUE; 3841 } else { 3842 b->free_a = PETSC_FALSE; 3843 b->free_ij = PETSC_FALSE; 3844 } 3845 3846 if (b->ipre && nnz != b->ipre && b->imax) { 3847 /* reserve user-requested sparsity */ 3848 ierr = PetscMemcpy(b->ipre,b->imax,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3849 } 3850 3851 3852 b->nz = 0; 3853 b->maxnz = nz; 3854 B->info.nz_unneeded = (double)b->maxnz; 3855 if (realalloc) { 3856 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3857 } 3858 B->was_assembled = PETSC_FALSE; 3859 B->assembled = PETSC_FALSE; 3860 PetscFunctionReturn(0); 3861 } 3862 3863 3864 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 3865 { 3866 Mat_SeqAIJ *a; 3867 PetscInt i; 3868 PetscErrorCode ierr; 3869 3870 PetscFunctionBegin; 3871 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3872 a = (Mat_SeqAIJ*)A->data; 3873 /* if no saved info, we error out */ 3874 if (!a->ipre) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"No saved preallocation info \n"); 3875 3876 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"); 3877 3878 ierr = PetscMemcpy(a->imax,a->ipre,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3879 ierr = PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3880 a->i[0] = 0; 3881 for (i=1; i<A->rmap->n+1; i++) { 3882 a->i[i] = a->i[i-1] + a->imax[i-1]; 3883 } 3884 A->preallocated = PETSC_TRUE; 3885 a->nz = 0; 3886 a->maxnz = a->i[A->rmap->n]; 3887 A->info.nz_unneeded = (double)a->maxnz; 3888 A->was_assembled = PETSC_FALSE; 3889 A->assembled = PETSC_FALSE; 3890 PetscFunctionReturn(0); 3891 } 3892 3893 /*@ 3894 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3895 3896 Input Parameters: 3897 + B - the matrix 3898 . i - the indices into j for the start of each row (starts with zero) 3899 . j - the column indices for each row (starts with zero) these must be sorted for each row 3900 - v - optional values in the matrix 3901 3902 Level: developer 3903 3904 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 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,NULL,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,NULL,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,NULL,PETSC_INT);CHKERRQ(ierr); 4458 4459 /* read in nonzero values */ 4460 ierr = PetscBinaryRead(fd,a->a,nz,NULL,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 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 4800 @*/ 4801 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 4802 { 4803 PetscErrorCode ierr,(*r)(Mat,MatType,MatReuse,Mat*); 4804 PetscBool sametype; 4805 4806 PetscFunctionBegin; 4807 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 4808 ierr = PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);CHKERRQ(ierr); 4809 if (sametype) PetscFunctionReturn(0); 4810 4811 ierr = PetscFunctionListFind(MatSeqAIJList,matype,&r);CHKERRQ(ierr); 4812 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 4813 ierr = (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);CHKERRQ(ierr); 4814 PetscFunctionReturn(0); 4815 } 4816 4817 4818 /*@C 4819 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 4820 4821 Not Collective 4822 4823 Input Parameters: 4824 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 4825 - function - routine to convert to subtype 4826 4827 Notes: 4828 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 4829 4830 4831 Then, your matrix can be chosen with the procedural interface at runtime via the option 4832 $ -mat_seqaij_type my_mat 4833 4834 Level: advanced 4835 4836 .seealso: MatSeqAIJRegisterAll() 4837 4838 4839 Level: advanced 4840 @*/ 4841 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 4842 { 4843 PetscErrorCode ierr; 4844 4845 PetscFunctionBegin; 4846 ierr = MatInitializePackage();CHKERRQ(ierr); 4847 ierr = PetscFunctionListAdd(&MatSeqAIJList,sname,function);CHKERRQ(ierr); 4848 PetscFunctionReturn(0); 4849 } 4850 4851 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 4852 4853 /*@C 4854 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 4855 4856 Not Collective 4857 4858 Level: advanced 4859 4860 Developers Note: CUSP and CUSPARSE do not yet support the MatConvert_SeqAIJ..() paradigm and thus cannot be registered here 4861 4862 .seealso: MatRegisterAll(), MatSeqAIJRegister() 4863 @*/ 4864 PetscErrorCode MatSeqAIJRegisterAll(void) 4865 { 4866 PetscErrorCode ierr; 4867 4868 PetscFunctionBegin; 4869 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 4870 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 4871 4872 ierr = MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4873 ierr = MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4874 ierr = MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL);CHKERRQ(ierr); 4875 #if defined(PETSC_HAVE_MKL_SPARSE) 4876 ierr = MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4877 #endif 4878 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 4879 ierr = MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr); 4880 #endif 4881 PetscFunctionReturn(0); 4882 } 4883 4884 /* 4885 Special version for direct calls from Fortran 4886 */ 4887 #include <petsc/private/fortranimpl.h> 4888 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4889 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4890 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4891 #define matsetvaluesseqaij_ matsetvaluesseqaij 4892 #endif 4893 4894 /* Change these macros so can be used in void function */ 4895 #undef CHKERRQ 4896 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4897 #undef SETERRQ2 4898 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4899 #undef SETERRQ3 4900 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4901 4902 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) 4903 { 4904 Mat A = *AA; 4905 PetscInt m = *mm, n = *nn; 4906 InsertMode is = *isis; 4907 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4908 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4909 PetscInt *imax,*ai,*ailen; 4910 PetscErrorCode ierr; 4911 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4912 MatScalar *ap,value,*aa; 4913 PetscBool ignorezeroentries = a->ignorezeroentries; 4914 PetscBool roworiented = a->roworiented; 4915 4916 PetscFunctionBegin; 4917 MatCheckPreallocated(A,1); 4918 imax = a->imax; 4919 ai = a->i; 4920 ailen = a->ilen; 4921 aj = a->j; 4922 aa = a->a; 4923 4924 for (k=0; k<m; k++) { /* loop over added rows */ 4925 row = im[k]; 4926 if (row < 0) continue; 4927 #if defined(PETSC_USE_DEBUG) 4928 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4929 #endif 4930 rp = aj + ai[row]; ap = aa + ai[row]; 4931 rmax = imax[row]; nrow = ailen[row]; 4932 low = 0; 4933 high = nrow; 4934 for (l=0; l<n; l++) { /* loop over added columns */ 4935 if (in[l] < 0) continue; 4936 #if defined(PETSC_USE_DEBUG) 4937 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4938 #endif 4939 col = in[l]; 4940 if (roworiented) value = v[l + k*n]; 4941 else value = v[k + l*m]; 4942 4943 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4944 4945 if (col <= lastcol) low = 0; 4946 else high = nrow; 4947 lastcol = col; 4948 while (high-low > 5) { 4949 t = (low+high)/2; 4950 if (rp[t] > col) high = t; 4951 else low = t; 4952 } 4953 for (i=low; i<high; i++) { 4954 if (rp[i] > col) break; 4955 if (rp[i] == col) { 4956 if (is == ADD_VALUES) ap[i] += value; 4957 else ap[i] = value; 4958 goto noinsert; 4959 } 4960 } 4961 if (value == 0.0 && ignorezeroentries) goto noinsert; 4962 if (nonew == 1) goto noinsert; 4963 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4964 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4965 N = nrow++ - 1; a->nz++; high++; 4966 /* shift up all the later entries in this row */ 4967 for (ii=N; ii>=i; ii--) { 4968 rp[ii+1] = rp[ii]; 4969 ap[ii+1] = ap[ii]; 4970 } 4971 rp[i] = col; 4972 ap[i] = value; 4973 A->nonzerostate++; 4974 noinsert:; 4975 low = i + 1; 4976 } 4977 ailen[row] = nrow; 4978 } 4979 PetscFunctionReturnVoid(); 4980 } 4981