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