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