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,"%%matrix complex general\n");CHKERRQ(ierr); 757 #else 758 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix 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 *x,*y; 1247 PetscErrorCode ierr; 1248 PetscInt m = A->rmap->n; 1249 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1250 MatScalar *v; 1251 PetscScalar alpha; 1252 PetscInt n,i,j,*idx,*ii,*ridx=NULL; 1253 Mat_CompressedRow cprow = a->compressedrow; 1254 PetscBool usecprow = cprow.use; 1255 #endif 1256 1257 PetscFunctionBegin; 1258 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 1259 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1260 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1261 1262 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1263 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 1264 #else 1265 if (usecprow) { 1266 m = cprow.nrows; 1267 ii = cprow.i; 1268 ridx = cprow.rindex; 1269 } else { 1270 ii = a->i; 1271 } 1272 for (i=0; i<m; i++) { 1273 idx = a->j + ii[i]; 1274 v = a->a + ii[i]; 1275 n = ii[i+1] - ii[i]; 1276 if (usecprow) { 1277 alpha = x[ridx[i]]; 1278 } else { 1279 alpha = x[i]; 1280 } 1281 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1282 } 1283 #endif 1284 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1285 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1286 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1287 PetscFunctionReturn(0); 1288 } 1289 1290 #undef __FUNCT__ 1291 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 1292 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1293 { 1294 PetscErrorCode ierr; 1295 1296 PetscFunctionBegin; 1297 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 1298 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 1299 PetscFunctionReturn(0); 1300 } 1301 1302 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1303 #if defined(PETSC_THREADCOMM_ACTIVE) 1304 PetscErrorCode MatMult_SeqAIJ_Kernel(PetscInt thread_id,Mat A,Vec xx,Vec yy) 1305 { 1306 PetscErrorCode ierr; 1307 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1308 PetscScalar *y; 1309 const PetscScalar *x; 1310 const MatScalar *aa; 1311 PetscInt *trstarts=A->rmap->trstarts; 1312 PetscInt n,start,end,i; 1313 const PetscInt *aj,*ai; 1314 PetscScalar sum; 1315 1316 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1317 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1318 start = trstarts[thread_id]; 1319 end = trstarts[thread_id+1]; 1320 aj = a->j; 1321 aa = a->a; 1322 ai = a->i; 1323 for (i=start; i<end; i++) { 1324 n = ai[i+1] - ai[i]; 1325 aj = a->j + ai[i]; 1326 aa = a->a + ai[i]; 1327 sum = 0.0; 1328 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1329 y[i] = sum; 1330 } 1331 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1332 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1333 return 0; 1334 } 1335 1336 #undef __FUNCT__ 1337 #define __FUNCT__ "MatMult_SeqAIJ" 1338 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1339 { 1340 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1341 PetscScalar *y; 1342 const PetscScalar *x; 1343 const MatScalar *aa; 1344 PetscErrorCode ierr; 1345 PetscInt m=A->rmap->n; 1346 const PetscInt *aj,*ii,*ridx=NULL; 1347 PetscInt n,i; 1348 PetscScalar sum; 1349 PetscBool usecprow=a->compressedrow.use; 1350 1351 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1352 #pragma disjoint(*x,*y,*aa) 1353 #endif 1354 1355 PetscFunctionBegin; 1356 aj = a->j; 1357 aa = a->a; 1358 ii = a->i; 1359 if (usecprow) { /* use compressed row format */ 1360 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1361 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1362 ierr = PetscMemzero(y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1363 m = a->compressedrow.nrows; 1364 ii = a->compressedrow.i; 1365 ridx = a->compressedrow.rindex; 1366 for (i=0; i<m; i++) { 1367 n = ii[i+1] - ii[i]; 1368 aj = a->j + ii[i]; 1369 aa = a->a + ii[i]; 1370 sum = 0.0; 1371 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1372 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1373 y[*ridx++] = sum; 1374 } 1375 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1376 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1377 } else { /* do not use compressed row format */ 1378 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1379 fortranmultaij_(&m,x,ii,aj,aa,y); 1380 #else 1381 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);CHKERRQ(ierr); 1382 #endif 1383 } 1384 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1385 PetscFunctionReturn(0); 1386 } 1387 #else 1388 #undef __FUNCT__ 1389 #define __FUNCT__ "MatMult_SeqAIJ" 1390 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1391 { 1392 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1393 PetscScalar *y; 1394 const PetscScalar *x; 1395 const MatScalar *aa; 1396 PetscErrorCode ierr; 1397 PetscInt m=A->rmap->n; 1398 const PetscInt *aj,*ii,*ridx=NULL; 1399 PetscInt n,i; 1400 PetscScalar sum; 1401 PetscBool usecprow=a->compressedrow.use; 1402 1403 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1404 #pragma disjoint(*x,*y,*aa) 1405 #endif 1406 1407 PetscFunctionBegin; 1408 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1409 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1410 aj = a->j; 1411 aa = a->a; 1412 ii = a->i; 1413 if (usecprow) { /* use compressed row format */ 1414 ierr = PetscMemzero(y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1415 m = a->compressedrow.nrows; 1416 ii = a->compressedrow.i; 1417 ridx = a->compressedrow.rindex; 1418 for (i=0; i<m; i++) { 1419 n = ii[i+1] - ii[i]; 1420 aj = a->j + ii[i]; 1421 aa = a->a + ii[i]; 1422 sum = 0.0; 1423 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1424 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1425 y[*ridx++] = sum; 1426 } 1427 } else { /* do not use compressed row format */ 1428 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1429 fortranmultaij_(&m,x,ii,aj,aa,y); 1430 #else 1431 #if defined(PETSC_THREADCOMM_ACTIVE) 1432 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);CHKERRQ(ierr); 1433 #else 1434 for (i=0; i<m; i++) { 1435 n = ii[i+1] - ii[i]; 1436 aj = a->j + ii[i]; 1437 aa = a->a + ii[i]; 1438 sum = 0.0; 1439 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1440 y[i] = sum; 1441 } 1442 #endif 1443 #endif 1444 } 1445 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1446 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1447 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1448 PetscFunctionReturn(0); 1449 } 1450 #endif 1451 1452 #undef __FUNCT__ 1453 #define __FUNCT__ "MatMultMax_SeqAIJ" 1454 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1455 { 1456 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1457 PetscScalar *y; 1458 const PetscScalar *x; 1459 const MatScalar *aa; 1460 PetscErrorCode ierr; 1461 PetscInt m=A->rmap->n; 1462 const PetscInt *aj,*ii,*ridx=NULL; 1463 PetscInt n,i,nonzerorow=0; 1464 PetscScalar sum; 1465 PetscBool usecprow=a->compressedrow.use; 1466 1467 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1468 #pragma disjoint(*x,*y,*aa) 1469 #endif 1470 1471 PetscFunctionBegin; 1472 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1473 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1474 aj = a->j; 1475 aa = a->a; 1476 ii = a->i; 1477 if (usecprow) { /* use compressed row format */ 1478 m = a->compressedrow.nrows; 1479 ii = a->compressedrow.i; 1480 ridx = a->compressedrow.rindex; 1481 for (i=0; i<m; i++) { 1482 n = ii[i+1] - ii[i]; 1483 aj = a->j + ii[i]; 1484 aa = a->a + ii[i]; 1485 sum = 0.0; 1486 nonzerorow += (n>0); 1487 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1488 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1489 y[*ridx++] = sum; 1490 } 1491 } else { /* do not use compressed row format */ 1492 for (i=0; i<m; i++) { 1493 n = ii[i+1] - ii[i]; 1494 aj = a->j + ii[i]; 1495 aa = a->a + ii[i]; 1496 sum = 0.0; 1497 nonzerorow += (n>0); 1498 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1499 y[i] = sum; 1500 } 1501 } 1502 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1503 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1504 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1505 PetscFunctionReturn(0); 1506 } 1507 1508 #undef __FUNCT__ 1509 #define __FUNCT__ "MatMultAddMax_SeqAIJ" 1510 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1511 { 1512 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1513 PetscScalar *y,*z; 1514 const PetscScalar *x; 1515 const MatScalar *aa; 1516 PetscErrorCode ierr; 1517 PetscInt m = A->rmap->n,*aj,*ii; 1518 PetscInt n,i,*ridx=NULL; 1519 PetscScalar sum; 1520 PetscBool usecprow=a->compressedrow.use; 1521 1522 PetscFunctionBegin; 1523 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1524 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1525 if (zz != yy) { 1526 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1527 } else { 1528 z = y; 1529 } 1530 1531 aj = a->j; 1532 aa = a->a; 1533 ii = a->i; 1534 if (usecprow) { /* use compressed row format */ 1535 if (zz != yy) { 1536 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1537 } 1538 m = a->compressedrow.nrows; 1539 ii = a->compressedrow.i; 1540 ridx = a->compressedrow.rindex; 1541 for (i=0; i<m; i++) { 1542 n = ii[i+1] - ii[i]; 1543 aj = a->j + ii[i]; 1544 aa = a->a + ii[i]; 1545 sum = y[*ridx]; 1546 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1547 z[*ridx++] = sum; 1548 } 1549 } else { /* do not use compressed row format */ 1550 for (i=0; i<m; i++) { 1551 n = ii[i+1] - ii[i]; 1552 aj = a->j + ii[i]; 1553 aa = a->a + ii[i]; 1554 sum = y[i]; 1555 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1556 z[i] = sum; 1557 } 1558 } 1559 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1560 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1561 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1562 if (zz != yy) { 1563 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1564 } 1565 PetscFunctionReturn(0); 1566 } 1567 1568 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1569 #undef __FUNCT__ 1570 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1571 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1572 { 1573 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1574 PetscScalar *y,*z; 1575 const PetscScalar *x; 1576 const MatScalar *aa; 1577 PetscErrorCode ierr; 1578 PetscInt m = A->rmap->n,*aj,*ii; 1579 PetscInt n,i,*ridx=NULL; 1580 PetscScalar sum; 1581 PetscBool usecprow=a->compressedrow.use; 1582 1583 PetscFunctionBegin; 1584 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1585 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1586 if (zz != yy) { 1587 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1588 } else { 1589 z = y; 1590 } 1591 1592 aj = a->j; 1593 aa = a->a; 1594 ii = a->i; 1595 if (usecprow) { /* use compressed row format */ 1596 if (zz != yy) { 1597 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1598 } 1599 m = a->compressedrow.nrows; 1600 ii = a->compressedrow.i; 1601 ridx = a->compressedrow.rindex; 1602 for (i=0; i<m; i++) { 1603 n = ii[i+1] - ii[i]; 1604 aj = a->j + ii[i]; 1605 aa = a->a + ii[i]; 1606 sum = y[*ridx]; 1607 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1608 z[*ridx++] = sum; 1609 } 1610 } else { /* do not use compressed row format */ 1611 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1612 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1613 #else 1614 for (i=0; i<m; i++) { 1615 n = ii[i+1] - ii[i]; 1616 aj = a->j + ii[i]; 1617 aa = a->a + ii[i]; 1618 sum = y[i]; 1619 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1620 z[i] = sum; 1621 } 1622 #endif 1623 } 1624 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1625 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1626 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1627 if (zz != yy) { 1628 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1629 } 1630 #if defined(PETSC_HAVE_CUSP) 1631 /* 1632 ierr = VecView(xx,0);CHKERRQ(ierr); 1633 ierr = VecView(zz,0);CHKERRQ(ierr); 1634 ierr = MatView(A,0);CHKERRQ(ierr); 1635 */ 1636 #endif 1637 PetscFunctionReturn(0); 1638 } 1639 1640 /* 1641 Adds diagonal pointers to sparse matrix structure. 1642 */ 1643 #undef __FUNCT__ 1644 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1645 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1646 { 1647 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1648 PetscErrorCode ierr; 1649 PetscInt i,j,m = A->rmap->n; 1650 1651 PetscFunctionBegin; 1652 if (!a->diag) { 1653 ierr = PetscMalloc1(m,&a->diag);CHKERRQ(ierr); 1654 ierr = PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));CHKERRQ(ierr); 1655 } 1656 for (i=0; i<A->rmap->n; i++) { 1657 a->diag[i] = a->i[i+1]; 1658 for (j=a->i[i]; j<a->i[i+1]; j++) { 1659 if (a->j[j] == i) { 1660 a->diag[i] = j; 1661 break; 1662 } 1663 } 1664 } 1665 PetscFunctionReturn(0); 1666 } 1667 1668 /* 1669 Checks for missing diagonals 1670 */ 1671 #undef __FUNCT__ 1672 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1673 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1674 { 1675 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1676 PetscInt *diag,*ii = a->i,i; 1677 1678 PetscFunctionBegin; 1679 *missing = PETSC_FALSE; 1680 if (A->rmap->n > 0 && !ii) { 1681 *missing = PETSC_TRUE; 1682 if (d) *d = 0; 1683 PetscInfo(A,"Matrix has no entries therefore is missing diagonal"); 1684 } else { 1685 diag = a->diag; 1686 for (i=0; i<A->rmap->n; i++) { 1687 if (diag[i] >= ii[i+1]) { 1688 *missing = PETSC_TRUE; 1689 if (d) *d = i; 1690 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1691 break; 1692 } 1693 } 1694 } 1695 PetscFunctionReturn(0); 1696 } 1697 1698 #undef __FUNCT__ 1699 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1700 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1701 { 1702 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1703 PetscErrorCode ierr; 1704 PetscInt i,*diag,m = A->rmap->n; 1705 MatScalar *v = a->a; 1706 PetscScalar *idiag,*mdiag; 1707 1708 PetscFunctionBegin; 1709 if (a->idiagvalid) PetscFunctionReturn(0); 1710 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1711 diag = a->diag; 1712 if (!a->idiag) { 1713 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1714 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1715 v = a->a; 1716 } 1717 mdiag = a->mdiag; 1718 idiag = a->idiag; 1719 1720 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1721 for (i=0; i<m; i++) { 1722 mdiag[i] = v[diag[i]]; 1723 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1724 idiag[i] = 1.0/v[diag[i]]; 1725 } 1726 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1727 } else { 1728 for (i=0; i<m; i++) { 1729 mdiag[i] = v[diag[i]]; 1730 idiag[i] = omega/(fshift + v[diag[i]]); 1731 } 1732 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1733 } 1734 a->idiagvalid = PETSC_TRUE; 1735 PetscFunctionReturn(0); 1736 } 1737 1738 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1739 #undef __FUNCT__ 1740 #define __FUNCT__ "MatSOR_SeqAIJ" 1741 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1742 { 1743 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1744 PetscScalar *x,d,sum,*t,scale; 1745 const MatScalar *v = a->a,*idiag=0,*mdiag; 1746 const PetscScalar *b, *bs,*xb, *ts; 1747 PetscErrorCode ierr; 1748 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1749 const PetscInt *idx,*diag; 1750 1751 PetscFunctionBegin; 1752 its = its*lits; 1753 1754 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1755 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1756 a->fshift = fshift; 1757 a->omega = omega; 1758 1759 diag = a->diag; 1760 t = a->ssor_work; 1761 idiag = a->idiag; 1762 mdiag = a->mdiag; 1763 1764 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1765 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1766 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1767 if (flag == SOR_APPLY_UPPER) { 1768 /* apply (U + D/omega) to the vector */ 1769 bs = b; 1770 for (i=0; i<m; i++) { 1771 d = fshift + mdiag[i]; 1772 n = a->i[i+1] - diag[i] - 1; 1773 idx = a->j + diag[i] + 1; 1774 v = a->a + diag[i] + 1; 1775 sum = b[i]*d/omega; 1776 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1777 x[i] = sum; 1778 } 1779 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1780 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1781 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1782 PetscFunctionReturn(0); 1783 } 1784 1785 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1786 else if (flag & SOR_EISENSTAT) { 1787 /* Let A = L + U + D; where L is lower trianglar, 1788 U is upper triangular, E = D/omega; This routine applies 1789 1790 (L + E)^{-1} A (U + E)^{-1} 1791 1792 to a vector efficiently using Eisenstat's trick. 1793 */ 1794 scale = (2.0/omega) - 1.0; 1795 1796 /* x = (E + U)^{-1} b */ 1797 for (i=m-1; i>=0; i--) { 1798 n = a->i[i+1] - diag[i] - 1; 1799 idx = a->j + diag[i] + 1; 1800 v = a->a + diag[i] + 1; 1801 sum = b[i]; 1802 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1803 x[i] = sum*idiag[i]; 1804 } 1805 1806 /* t = b - (2*E - D)x */ 1807 v = a->a; 1808 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1809 1810 /* t = (E + L)^{-1}t */ 1811 ts = t; 1812 diag = a->diag; 1813 for (i=0; i<m; i++) { 1814 n = diag[i] - a->i[i]; 1815 idx = a->j + a->i[i]; 1816 v = a->a + a->i[i]; 1817 sum = t[i]; 1818 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1819 t[i] = sum*idiag[i]; 1820 /* x = x + t */ 1821 x[i] += t[i]; 1822 } 1823 1824 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1825 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1826 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1827 PetscFunctionReturn(0); 1828 } 1829 if (flag & SOR_ZERO_INITIAL_GUESS) { 1830 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1831 for (i=0; i<m; i++) { 1832 n = diag[i] - a->i[i]; 1833 idx = a->j + a->i[i]; 1834 v = a->a + a->i[i]; 1835 sum = b[i]; 1836 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1837 t[i] = sum; 1838 x[i] = sum*idiag[i]; 1839 } 1840 xb = t; 1841 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1842 } else xb = b; 1843 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1844 for (i=m-1; i>=0; i--) { 1845 n = a->i[i+1] - diag[i] - 1; 1846 idx = a->j + diag[i] + 1; 1847 v = a->a + diag[i] + 1; 1848 sum = xb[i]; 1849 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1850 if (xb == b) { 1851 x[i] = sum*idiag[i]; 1852 } else { 1853 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1854 } 1855 } 1856 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1857 } 1858 its--; 1859 } 1860 while (its--) { 1861 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1862 for (i=0; i<m; i++) { 1863 /* lower */ 1864 n = diag[i] - a->i[i]; 1865 idx = a->j + a->i[i]; 1866 v = a->a + a->i[i]; 1867 sum = b[i]; 1868 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1869 t[i] = sum; /* save application of the lower-triangular part */ 1870 /* upper */ 1871 n = a->i[i+1] - diag[i] - 1; 1872 idx = a->j + diag[i] + 1; 1873 v = a->a + diag[i] + 1; 1874 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1875 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1876 } 1877 xb = t; 1878 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1879 } else xb = b; 1880 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1881 for (i=m-1; i>=0; i--) { 1882 sum = xb[i]; 1883 if (xb == b) { 1884 /* whole matrix (no checkpointing available) */ 1885 n = a->i[i+1] - a->i[i]; 1886 idx = a->j + a->i[i]; 1887 v = a->a + a->i[i]; 1888 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1889 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1890 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1891 n = a->i[i+1] - diag[i] - 1; 1892 idx = a->j + diag[i] + 1; 1893 v = a->a + diag[i] + 1; 1894 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1895 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1896 } 1897 } 1898 if (xb == b) { 1899 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1900 } else { 1901 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1902 } 1903 } 1904 } 1905 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1906 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1907 PetscFunctionReturn(0); 1908 } 1909 1910 1911 #undef __FUNCT__ 1912 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1913 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1914 { 1915 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1916 1917 PetscFunctionBegin; 1918 info->block_size = 1.0; 1919 info->nz_allocated = (double)a->maxnz; 1920 info->nz_used = (double)a->nz; 1921 info->nz_unneeded = (double)(a->maxnz - a->nz); 1922 info->assemblies = (double)A->num_ass; 1923 info->mallocs = (double)A->info.mallocs; 1924 info->memory = ((PetscObject)A)->mem; 1925 if (A->factortype) { 1926 info->fill_ratio_given = A->info.fill_ratio_given; 1927 info->fill_ratio_needed = A->info.fill_ratio_needed; 1928 info->factor_mallocs = A->info.factor_mallocs; 1929 } else { 1930 info->fill_ratio_given = 0; 1931 info->fill_ratio_needed = 0; 1932 info->factor_mallocs = 0; 1933 } 1934 PetscFunctionReturn(0); 1935 } 1936 1937 #undef __FUNCT__ 1938 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1939 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1940 { 1941 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1942 PetscInt i,m = A->rmap->n - 1,d = 0; 1943 PetscErrorCode ierr; 1944 const PetscScalar *xx; 1945 PetscScalar *bb; 1946 PetscBool missing; 1947 1948 PetscFunctionBegin; 1949 if (x && b) { 1950 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1951 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1952 for (i=0; i<N; i++) { 1953 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1954 bb[rows[i]] = diag*xx[rows[i]]; 1955 } 1956 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1957 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1958 } 1959 1960 if (a->keepnonzeropattern) { 1961 for (i=0; i<N; i++) { 1962 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1963 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1964 } 1965 if (diag != 0.0) { 1966 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1967 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1968 for (i=0; i<N; i++) { 1969 a->a[a->diag[rows[i]]] = diag; 1970 } 1971 } 1972 } else { 1973 if (diag != 0.0) { 1974 for (i=0; i<N; i++) { 1975 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1976 if (a->ilen[rows[i]] > 0) { 1977 a->ilen[rows[i]] = 1; 1978 a->a[a->i[rows[i]]] = diag; 1979 a->j[a->i[rows[i]]] = rows[i]; 1980 } else { /* in case row was completely empty */ 1981 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1982 } 1983 } 1984 } else { 1985 for (i=0; i<N; i++) { 1986 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1987 a->ilen[rows[i]] = 0; 1988 } 1989 } 1990 A->nonzerostate++; 1991 } 1992 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1993 PetscFunctionReturn(0); 1994 } 1995 1996 #undef __FUNCT__ 1997 #define __FUNCT__ "MatZeroRowsColumns_SeqAIJ" 1998 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1999 { 2000 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2001 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2002 PetscErrorCode ierr; 2003 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2004 const PetscScalar *xx; 2005 PetscScalar *bb; 2006 2007 PetscFunctionBegin; 2008 if (x && b) { 2009 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2010 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 2011 vecs = PETSC_TRUE; 2012 } 2013 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 2014 for (i=0; i<N; i++) { 2015 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 2016 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 2017 2018 zeroed[rows[i]] = PETSC_TRUE; 2019 } 2020 for (i=0; i<A->rmap->n; i++) { 2021 if (!zeroed[i]) { 2022 for (j=a->i[i]; j<a->i[i+1]; j++) { 2023 if (zeroed[a->j[j]]) { 2024 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 2025 a->a[j] = 0.0; 2026 } 2027 } 2028 } else if (vecs) bb[i] = diag*xx[i]; 2029 } 2030 if (x && b) { 2031 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2032 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 2033 } 2034 ierr = PetscFree(zeroed);CHKERRQ(ierr); 2035 if (diag != 0.0) { 2036 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 2037 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 2038 for (i=0; i<N; i++) { 2039 a->a[a->diag[rows[i]]] = diag; 2040 } 2041 } 2042 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2043 PetscFunctionReturn(0); 2044 } 2045 2046 #undef __FUNCT__ 2047 #define __FUNCT__ "MatGetRow_SeqAIJ" 2048 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2049 { 2050 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2051 PetscInt *itmp; 2052 2053 PetscFunctionBegin; 2054 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 2055 2056 *nz = a->i[row+1] - a->i[row]; 2057 if (v) *v = a->a + a->i[row]; 2058 if (idx) { 2059 itmp = a->j + a->i[row]; 2060 if (*nz) *idx = itmp; 2061 else *idx = 0; 2062 } 2063 PetscFunctionReturn(0); 2064 } 2065 2066 /* remove this function? */ 2067 #undef __FUNCT__ 2068 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 2069 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2070 { 2071 PetscFunctionBegin; 2072 PetscFunctionReturn(0); 2073 } 2074 2075 #undef __FUNCT__ 2076 #define __FUNCT__ "MatNorm_SeqAIJ" 2077 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2078 { 2079 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2080 MatScalar *v = a->a; 2081 PetscReal sum = 0.0; 2082 PetscErrorCode ierr; 2083 PetscInt i,j; 2084 2085 PetscFunctionBegin; 2086 if (type == NORM_FROBENIUS) { 2087 for (i=0; i<a->nz; i++) { 2088 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2089 } 2090 *nrm = PetscSqrtReal(sum); 2091 } else if (type == NORM_1) { 2092 PetscReal *tmp; 2093 PetscInt *jj = a->j; 2094 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 2095 *nrm = 0.0; 2096 for (j=0; j<a->nz; j++) { 2097 tmp[*jj++] += PetscAbsScalar(*v); v++; 2098 } 2099 for (j=0; j<A->cmap->n; j++) { 2100 if (tmp[j] > *nrm) *nrm = tmp[j]; 2101 } 2102 ierr = PetscFree(tmp);CHKERRQ(ierr); 2103 } else if (type == NORM_INFINITY) { 2104 *nrm = 0.0; 2105 for (j=0; j<A->rmap->n; j++) { 2106 v = a->a + a->i[j]; 2107 sum = 0.0; 2108 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2109 sum += PetscAbsScalar(*v); v++; 2110 } 2111 if (sum > *nrm) *nrm = sum; 2112 } 2113 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2114 PetscFunctionReturn(0); 2115 } 2116 2117 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2118 #undef __FUNCT__ 2119 #define __FUNCT__ "MatTransposeSymbolic_SeqAIJ" 2120 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2121 { 2122 PetscErrorCode ierr; 2123 PetscInt i,j,anzj; 2124 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2125 PetscInt an=A->cmap->N,am=A->rmap->N; 2126 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2127 2128 PetscFunctionBegin; 2129 /* Allocate space for symbolic transpose info and work array */ 2130 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 2131 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 2132 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 2133 2134 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2135 /* Note: offset by 1 for fast conversion into csr format. */ 2136 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2137 /* Form ati for csr format of A^T. */ 2138 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2139 2140 /* Copy ati into atfill so we have locations of the next free space in atj */ 2141 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2142 2143 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2144 for (i=0;i<am;i++) { 2145 anzj = ai[i+1] - ai[i]; 2146 for (j=0;j<anzj;j++) { 2147 atj[atfill[*aj]] = i; 2148 atfill[*aj++] += 1; 2149 } 2150 } 2151 2152 /* Clean up temporary space and complete requests. */ 2153 ierr = PetscFree(atfill);CHKERRQ(ierr); 2154 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2155 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2156 2157 b = (Mat_SeqAIJ*)((*B)->data); 2158 b->free_a = PETSC_FALSE; 2159 b->free_ij = PETSC_TRUE; 2160 b->nonew = 0; 2161 PetscFunctionReturn(0); 2162 } 2163 2164 #undef __FUNCT__ 2165 #define __FUNCT__ "MatTranspose_SeqAIJ" 2166 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2167 { 2168 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2169 Mat C; 2170 PetscErrorCode ierr; 2171 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2172 MatScalar *array = a->a; 2173 2174 PetscFunctionBegin; 2175 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"); 2176 2177 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 2178 ierr = PetscCalloc1(1+A->cmap->n,&col);CHKERRQ(ierr); 2179 2180 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2181 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2182 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2183 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2184 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2185 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2186 ierr = PetscFree(col);CHKERRQ(ierr); 2187 } else { 2188 C = *B; 2189 } 2190 2191 for (i=0; i<m; i++) { 2192 len = ai[i+1]-ai[i]; 2193 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2194 array += len; 2195 aj += len; 2196 } 2197 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2198 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2199 2200 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 2201 *B = C; 2202 } else { 2203 ierr = MatHeaderMerge(A,C);CHKERRQ(ierr); 2204 } 2205 PetscFunctionReturn(0); 2206 } 2207 2208 #undef __FUNCT__ 2209 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 2210 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2211 { 2212 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2213 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2214 MatScalar *va,*vb; 2215 PetscErrorCode ierr; 2216 PetscInt ma,na,mb,nb, i; 2217 2218 PetscFunctionBegin; 2219 bij = (Mat_SeqAIJ*) B->data; 2220 2221 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2222 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2223 if (ma!=nb || na!=mb) { 2224 *f = PETSC_FALSE; 2225 PetscFunctionReturn(0); 2226 } 2227 aii = aij->i; bii = bij->i; 2228 adx = aij->j; bdx = bij->j; 2229 va = aij->a; vb = bij->a; 2230 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2231 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2232 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2233 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2234 2235 *f = PETSC_TRUE; 2236 for (i=0; i<ma; i++) { 2237 while (aptr[i]<aii[i+1]) { 2238 PetscInt idc,idr; 2239 PetscScalar vc,vr; 2240 /* column/row index/value */ 2241 idc = adx[aptr[i]]; 2242 idr = bdx[bptr[idc]]; 2243 vc = va[aptr[i]]; 2244 vr = vb[bptr[idc]]; 2245 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2246 *f = PETSC_FALSE; 2247 goto done; 2248 } else { 2249 aptr[i]++; 2250 if (B || i!=idc) bptr[idc]++; 2251 } 2252 } 2253 } 2254 done: 2255 ierr = PetscFree(aptr);CHKERRQ(ierr); 2256 ierr = PetscFree(bptr);CHKERRQ(ierr); 2257 PetscFunctionReturn(0); 2258 } 2259 2260 #undef __FUNCT__ 2261 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 2262 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2263 { 2264 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2265 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2266 MatScalar *va,*vb; 2267 PetscErrorCode ierr; 2268 PetscInt ma,na,mb,nb, i; 2269 2270 PetscFunctionBegin; 2271 bij = (Mat_SeqAIJ*) B->data; 2272 2273 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2274 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2275 if (ma!=nb || na!=mb) { 2276 *f = PETSC_FALSE; 2277 PetscFunctionReturn(0); 2278 } 2279 aii = aij->i; bii = bij->i; 2280 adx = aij->j; bdx = bij->j; 2281 va = aij->a; vb = bij->a; 2282 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2283 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2284 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2285 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2286 2287 *f = PETSC_TRUE; 2288 for (i=0; i<ma; i++) { 2289 while (aptr[i]<aii[i+1]) { 2290 PetscInt idc,idr; 2291 PetscScalar vc,vr; 2292 /* column/row index/value */ 2293 idc = adx[aptr[i]]; 2294 idr = bdx[bptr[idc]]; 2295 vc = va[aptr[i]]; 2296 vr = vb[bptr[idc]]; 2297 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2298 *f = PETSC_FALSE; 2299 goto done; 2300 } else { 2301 aptr[i]++; 2302 if (B || i!=idc) bptr[idc]++; 2303 } 2304 } 2305 } 2306 done: 2307 ierr = PetscFree(aptr);CHKERRQ(ierr); 2308 ierr = PetscFree(bptr);CHKERRQ(ierr); 2309 PetscFunctionReturn(0); 2310 } 2311 2312 #undef __FUNCT__ 2313 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 2314 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2315 { 2316 PetscErrorCode ierr; 2317 2318 PetscFunctionBegin; 2319 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2320 PetscFunctionReturn(0); 2321 } 2322 2323 #undef __FUNCT__ 2324 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 2325 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2326 { 2327 PetscErrorCode ierr; 2328 2329 PetscFunctionBegin; 2330 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2331 PetscFunctionReturn(0); 2332 } 2333 2334 #undef __FUNCT__ 2335 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 2336 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2337 { 2338 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2339 PetscScalar *l,*r,x; 2340 MatScalar *v; 2341 PetscErrorCode ierr; 2342 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 2343 2344 PetscFunctionBegin; 2345 if (ll) { 2346 /* The local size is used so that VecMPI can be passed to this routine 2347 by MatDiagonalScale_MPIAIJ */ 2348 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2349 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2350 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 2351 v = a->a; 2352 for (i=0; i<m; i++) { 2353 x = l[i]; 2354 M = a->i[i+1] - a->i[i]; 2355 for (j=0; j<M; j++) (*v++) *= x; 2356 } 2357 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 2358 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2359 } 2360 if (rr) { 2361 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2362 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2363 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 2364 v = a->a; jj = a->j; 2365 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2366 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 2367 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2368 } 2369 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2370 PetscFunctionReturn(0); 2371 } 2372 2373 #undef __FUNCT__ 2374 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 2375 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2376 { 2377 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2378 PetscErrorCode ierr; 2379 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2380 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2381 const PetscInt *irow,*icol; 2382 PetscInt nrows,ncols; 2383 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2384 MatScalar *a_new,*mat_a; 2385 Mat C; 2386 PetscBool stride,sorted; 2387 2388 PetscFunctionBegin; 2389 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 2390 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 2391 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 2392 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 2393 2394 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2395 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2396 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2397 2398 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2399 if (stride) { 2400 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2401 } else { 2402 first = 0; 2403 step = 0; 2404 } 2405 if (stride && step == 1) { 2406 /* special case of contiguous rows */ 2407 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2408 /* loop over new rows determining lens and starting points */ 2409 for (i=0; i<nrows; i++) { 2410 kstart = ai[irow[i]]; 2411 kend = kstart + ailen[irow[i]]; 2412 starts[i] = kstart; 2413 for (k=kstart; k<kend; k++) { 2414 if (aj[k] >= first) { 2415 starts[i] = k; 2416 break; 2417 } 2418 } 2419 sum = 0; 2420 while (k < kend) { 2421 if (aj[k++] >= first+ncols) break; 2422 sum++; 2423 } 2424 lens[i] = sum; 2425 } 2426 /* create submatrix */ 2427 if (scall == MAT_REUSE_MATRIX) { 2428 PetscInt n_cols,n_rows; 2429 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2430 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2431 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2432 C = *B; 2433 } else { 2434 PetscInt rbs,cbs; 2435 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2436 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2437 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2438 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2439 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2440 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2441 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2442 } 2443 c = (Mat_SeqAIJ*)C->data; 2444 2445 /* loop over rows inserting into submatrix */ 2446 a_new = c->a; 2447 j_new = c->j; 2448 i_new = c->i; 2449 2450 for (i=0; i<nrows; i++) { 2451 ii = starts[i]; 2452 lensi = lens[i]; 2453 for (k=0; k<lensi; k++) { 2454 *j_new++ = aj[ii+k] - first; 2455 } 2456 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2457 a_new += lensi; 2458 i_new[i+1] = i_new[i] + lensi; 2459 c->ilen[i] = lensi; 2460 } 2461 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2462 } else { 2463 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2464 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2465 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2466 for (i=0; i<ncols; i++) { 2467 #if defined(PETSC_USE_DEBUG) 2468 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); 2469 #endif 2470 smap[icol[i]] = i+1; 2471 } 2472 2473 /* determine lens of each row */ 2474 for (i=0; i<nrows; i++) { 2475 kstart = ai[irow[i]]; 2476 kend = kstart + a->ilen[irow[i]]; 2477 lens[i] = 0; 2478 for (k=kstart; k<kend; k++) { 2479 if (smap[aj[k]]) { 2480 lens[i]++; 2481 } 2482 } 2483 } 2484 /* Create and fill new matrix */ 2485 if (scall == MAT_REUSE_MATRIX) { 2486 PetscBool equal; 2487 2488 c = (Mat_SeqAIJ*)((*B)->data); 2489 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2490 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2491 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2492 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2493 C = *B; 2494 } else { 2495 PetscInt rbs,cbs; 2496 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2497 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2498 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2499 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2500 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2501 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2502 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2503 } 2504 c = (Mat_SeqAIJ*)(C->data); 2505 for (i=0; i<nrows; i++) { 2506 row = irow[i]; 2507 kstart = ai[row]; 2508 kend = kstart + a->ilen[row]; 2509 mat_i = c->i[i]; 2510 mat_j = c->j + mat_i; 2511 mat_a = c->a + mat_i; 2512 mat_ilen = c->ilen + i; 2513 for (k=kstart; k<kend; k++) { 2514 if ((tcol=smap[a->j[k]])) { 2515 *mat_j++ = tcol - 1; 2516 *mat_a++ = a->a[k]; 2517 (*mat_ilen)++; 2518 2519 } 2520 } 2521 } 2522 /* Free work space */ 2523 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2524 ierr = PetscFree(smap);CHKERRQ(ierr); 2525 ierr = PetscFree(lens);CHKERRQ(ierr); 2526 } 2527 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2528 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2529 2530 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2531 *B = C; 2532 PetscFunctionReturn(0); 2533 } 2534 2535 #undef __FUNCT__ 2536 #define __FUNCT__ "MatGetMultiProcBlock_SeqAIJ" 2537 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2538 { 2539 PetscErrorCode ierr; 2540 Mat B; 2541 2542 PetscFunctionBegin; 2543 if (scall == MAT_INITIAL_MATRIX) { 2544 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2545 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2546 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2547 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2548 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2549 *subMat = B; 2550 } else { 2551 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2552 } 2553 PetscFunctionReturn(0); 2554 } 2555 2556 #undef __FUNCT__ 2557 #define __FUNCT__ "MatILUFactor_SeqAIJ" 2558 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2559 { 2560 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2561 PetscErrorCode ierr; 2562 Mat outA; 2563 PetscBool row_identity,col_identity; 2564 2565 PetscFunctionBegin; 2566 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2567 2568 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2569 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2570 2571 outA = inA; 2572 outA->factortype = MAT_FACTOR_LU; 2573 2574 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2575 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2576 2577 a->row = row; 2578 2579 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2580 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2581 2582 a->col = col; 2583 2584 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2585 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2586 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2587 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2588 2589 if (!a->solve_work) { /* this matrix may have been factored before */ 2590 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2591 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2592 } 2593 2594 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2595 if (row_identity && col_identity) { 2596 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2597 } else { 2598 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2599 } 2600 PetscFunctionReturn(0); 2601 } 2602 2603 #undef __FUNCT__ 2604 #define __FUNCT__ "MatScale_SeqAIJ" 2605 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2606 { 2607 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2608 PetscScalar oalpha = alpha; 2609 PetscErrorCode ierr; 2610 PetscBLASInt one = 1,bnz; 2611 2612 PetscFunctionBegin; 2613 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2614 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2615 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2616 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2617 PetscFunctionReturn(0); 2618 } 2619 2620 #undef __FUNCT__ 2621 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 2622 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2623 { 2624 PetscErrorCode ierr; 2625 PetscInt i; 2626 2627 PetscFunctionBegin; 2628 if (scall == MAT_INITIAL_MATRIX) { 2629 ierr = PetscMalloc1(n+1,B);CHKERRQ(ierr); 2630 } 2631 2632 for (i=0; i<n; i++) { 2633 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2634 } 2635 PetscFunctionReturn(0); 2636 } 2637 2638 #undef __FUNCT__ 2639 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 2640 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2641 { 2642 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2643 PetscErrorCode ierr; 2644 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2645 const PetscInt *idx; 2646 PetscInt start,end,*ai,*aj; 2647 PetscBT table; 2648 2649 PetscFunctionBegin; 2650 m = A->rmap->n; 2651 ai = a->i; 2652 aj = a->j; 2653 2654 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2655 2656 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2657 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2658 2659 for (i=0; i<is_max; i++) { 2660 /* Initialize the two local arrays */ 2661 isz = 0; 2662 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2663 2664 /* Extract the indices, assume there can be duplicate entries */ 2665 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2666 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2667 2668 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2669 for (j=0; j<n; ++j) { 2670 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2671 } 2672 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2673 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2674 2675 k = 0; 2676 for (j=0; j<ov; j++) { /* for each overlap */ 2677 n = isz; 2678 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2679 row = nidx[k]; 2680 start = ai[row]; 2681 end = ai[row+1]; 2682 for (l = start; l<end; l++) { 2683 val = aj[l]; 2684 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2685 } 2686 } 2687 } 2688 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2689 } 2690 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2691 ierr = PetscFree(nidx);CHKERRQ(ierr); 2692 PetscFunctionReturn(0); 2693 } 2694 2695 /* -------------------------------------------------------------- */ 2696 #undef __FUNCT__ 2697 #define __FUNCT__ "MatPermute_SeqAIJ" 2698 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2699 { 2700 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2701 PetscErrorCode ierr; 2702 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2703 const PetscInt *row,*col; 2704 PetscInt *cnew,j,*lens; 2705 IS icolp,irowp; 2706 PetscInt *cwork = NULL; 2707 PetscScalar *vwork = NULL; 2708 2709 PetscFunctionBegin; 2710 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2711 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2712 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2713 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2714 2715 /* determine lengths of permuted rows */ 2716 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2717 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2718 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2719 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2720 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2721 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2722 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2723 ierr = PetscFree(lens);CHKERRQ(ierr); 2724 2725 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2726 for (i=0; i<m; i++) { 2727 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2728 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2729 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2730 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2731 } 2732 ierr = PetscFree(cnew);CHKERRQ(ierr); 2733 2734 (*B)->assembled = PETSC_FALSE; 2735 2736 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2737 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2738 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2739 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2740 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2741 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2742 PetscFunctionReturn(0); 2743 } 2744 2745 #undef __FUNCT__ 2746 #define __FUNCT__ "MatCopy_SeqAIJ" 2747 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2748 { 2749 PetscErrorCode ierr; 2750 2751 PetscFunctionBegin; 2752 /* If the two matrices have the same copy implementation, use fast copy. */ 2753 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2754 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2755 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2756 2757 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"); 2758 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2759 } else { 2760 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2761 } 2762 PetscFunctionReturn(0); 2763 } 2764 2765 #undef __FUNCT__ 2766 #define __FUNCT__ "MatSetUp_SeqAIJ" 2767 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2768 { 2769 PetscErrorCode ierr; 2770 2771 PetscFunctionBegin; 2772 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2773 PetscFunctionReturn(0); 2774 } 2775 2776 #undef __FUNCT__ 2777 #define __FUNCT__ "MatSeqAIJGetArray_SeqAIJ" 2778 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2779 { 2780 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2781 2782 PetscFunctionBegin; 2783 *array = a->a; 2784 PetscFunctionReturn(0); 2785 } 2786 2787 #undef __FUNCT__ 2788 #define __FUNCT__ "MatSeqAIJRestoreArray_SeqAIJ" 2789 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2790 { 2791 PetscFunctionBegin; 2792 PetscFunctionReturn(0); 2793 } 2794 2795 /* 2796 Computes the number of nonzeros per row needed for preallocation when X and Y 2797 have different nonzero structure. 2798 */ 2799 #undef __FUNCT__ 2800 #define __FUNCT__ "MatAXPYGetPreallocation_SeqX_private" 2801 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2802 { 2803 PetscInt i,j,k,nzx,nzy; 2804 2805 PetscFunctionBegin; 2806 /* Set the number of nonzeros in the new matrix */ 2807 for (i=0; i<m; i++) { 2808 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2809 nzx = xi[i+1] - xi[i]; 2810 nzy = yi[i+1] - yi[i]; 2811 nnz[i] = 0; 2812 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2813 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2814 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2815 nnz[i]++; 2816 } 2817 for (; k<nzy; k++) nnz[i]++; 2818 } 2819 PetscFunctionReturn(0); 2820 } 2821 2822 #undef __FUNCT__ 2823 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2824 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2825 { 2826 PetscInt m = Y->rmap->N; 2827 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2828 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2829 PetscErrorCode ierr; 2830 2831 PetscFunctionBegin; 2832 /* Set the number of nonzeros in the new matrix */ 2833 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2834 PetscFunctionReturn(0); 2835 } 2836 2837 #undef __FUNCT__ 2838 #define __FUNCT__ "MatAXPY_SeqAIJ" 2839 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2840 { 2841 PetscErrorCode ierr; 2842 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2843 PetscBLASInt one=1,bnz; 2844 2845 PetscFunctionBegin; 2846 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2847 if (str == SAME_NONZERO_PATTERN) { 2848 PetscScalar alpha = a; 2849 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2850 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2851 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2852 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2853 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2854 } else { 2855 Mat B; 2856 PetscInt *nnz; 2857 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2858 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2859 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2860 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2861 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2862 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2863 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2864 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2865 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2866 ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr); 2867 ierr = PetscFree(nnz);CHKERRQ(ierr); 2868 } 2869 PetscFunctionReturn(0); 2870 } 2871 2872 #undef __FUNCT__ 2873 #define __FUNCT__ "MatConjugate_SeqAIJ" 2874 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2875 { 2876 #if defined(PETSC_USE_COMPLEX) 2877 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2878 PetscInt i,nz; 2879 PetscScalar *a; 2880 2881 PetscFunctionBegin; 2882 nz = aij->nz; 2883 a = aij->a; 2884 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2885 #else 2886 PetscFunctionBegin; 2887 #endif 2888 PetscFunctionReturn(0); 2889 } 2890 2891 #undef __FUNCT__ 2892 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2893 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2894 { 2895 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2896 PetscErrorCode ierr; 2897 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2898 PetscReal atmp; 2899 PetscScalar *x; 2900 MatScalar *aa; 2901 2902 PetscFunctionBegin; 2903 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2904 aa = a->a; 2905 ai = a->i; 2906 aj = a->j; 2907 2908 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2909 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2910 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2911 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2912 for (i=0; i<m; i++) { 2913 ncols = ai[1] - ai[0]; ai++; 2914 x[i] = 0.0; 2915 for (j=0; j<ncols; j++) { 2916 atmp = PetscAbsScalar(*aa); 2917 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2918 aa++; aj++; 2919 } 2920 } 2921 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2922 PetscFunctionReturn(0); 2923 } 2924 2925 #undef __FUNCT__ 2926 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2927 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2928 { 2929 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2930 PetscErrorCode ierr; 2931 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2932 PetscScalar *x; 2933 MatScalar *aa; 2934 2935 PetscFunctionBegin; 2936 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2937 aa = a->a; 2938 ai = a->i; 2939 aj = a->j; 2940 2941 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2942 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2943 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2944 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2945 for (i=0; i<m; i++) { 2946 ncols = ai[1] - ai[0]; ai++; 2947 if (ncols == A->cmap->n) { /* row is dense */ 2948 x[i] = *aa; if (idx) idx[i] = 0; 2949 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2950 x[i] = 0.0; 2951 if (idx) { 2952 idx[i] = 0; /* in case ncols is zero */ 2953 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2954 if (aj[j] > j) { 2955 idx[i] = j; 2956 break; 2957 } 2958 } 2959 } 2960 } 2961 for (j=0; j<ncols; j++) { 2962 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2963 aa++; aj++; 2964 } 2965 } 2966 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2967 PetscFunctionReturn(0); 2968 } 2969 2970 #undef __FUNCT__ 2971 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2972 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2973 { 2974 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2975 PetscErrorCode ierr; 2976 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2977 PetscReal atmp; 2978 PetscScalar *x; 2979 MatScalar *aa; 2980 2981 PetscFunctionBegin; 2982 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2983 aa = a->a; 2984 ai = a->i; 2985 aj = a->j; 2986 2987 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2988 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2989 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2990 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); 2991 for (i=0; i<m; i++) { 2992 ncols = ai[1] - ai[0]; ai++; 2993 if (ncols) { 2994 /* Get first nonzero */ 2995 for (j = 0; j < ncols; j++) { 2996 atmp = PetscAbsScalar(aa[j]); 2997 if (atmp > 1.0e-12) { 2998 x[i] = atmp; 2999 if (idx) idx[i] = aj[j]; 3000 break; 3001 } 3002 } 3003 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 3004 } else { 3005 x[i] = 0.0; if (idx) idx[i] = 0; 3006 } 3007 for (j = 0; j < ncols; j++) { 3008 atmp = PetscAbsScalar(*aa); 3009 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3010 aa++; aj++; 3011 } 3012 } 3013 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3014 PetscFunctionReturn(0); 3015 } 3016 3017 #undef __FUNCT__ 3018 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 3019 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3020 { 3021 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3022 PetscErrorCode ierr; 3023 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3024 PetscScalar *x; 3025 MatScalar *aa; 3026 3027 PetscFunctionBegin; 3028 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3029 aa = a->a; 3030 ai = a->i; 3031 aj = a->j; 3032 3033 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3034 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3035 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3036 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3037 for (i=0; i<m; i++) { 3038 ncols = ai[1] - ai[0]; ai++; 3039 if (ncols == A->cmap->n) { /* row is dense */ 3040 x[i] = *aa; if (idx) idx[i] = 0; 3041 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3042 x[i] = 0.0; 3043 if (idx) { /* find first implicit 0.0 in the row */ 3044 idx[i] = 0; /* in case ncols is zero */ 3045 for (j=0; j<ncols; j++) { 3046 if (aj[j] > j) { 3047 idx[i] = j; 3048 break; 3049 } 3050 } 3051 } 3052 } 3053 for (j=0; j<ncols; j++) { 3054 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3055 aa++; aj++; 3056 } 3057 } 3058 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3059 PetscFunctionReturn(0); 3060 } 3061 3062 #include <petscblaslapack.h> 3063 #include <petsc-private/kernels/blockinvert.h> 3064 3065 #undef __FUNCT__ 3066 #define __FUNCT__ "MatInvertBlockDiagonal_SeqAIJ" 3067 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3068 { 3069 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3070 PetscErrorCode ierr; 3071 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3072 MatScalar *diag,work[25],*v_work; 3073 PetscReal shift = 0.0; 3074 3075 PetscFunctionBegin; 3076 if (a->ibdiagvalid) { 3077 if (values) *values = a->ibdiag; 3078 PetscFunctionReturn(0); 3079 } 3080 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 3081 if (!a->ibdiag) { 3082 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 3083 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 3084 } 3085 diag = a->ibdiag; 3086 if (values) *values = a->ibdiag; 3087 /* factor and invert each block */ 3088 switch (bs) { 3089 case 1: 3090 for (i=0; i<mbs; i++) { 3091 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3092 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3093 } 3094 break; 3095 case 2: 3096 for (i=0; i<mbs; i++) { 3097 ij[0] = 2*i; ij[1] = 2*i + 1; 3098 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3099 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift);CHKERRQ(ierr); 3100 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3101 diag += 4; 3102 } 3103 break; 3104 case 3: 3105 for (i=0; i<mbs; i++) { 3106 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3107 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3108 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift);CHKERRQ(ierr); 3109 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3110 diag += 9; 3111 } 3112 break; 3113 case 4: 3114 for (i=0; i<mbs; i++) { 3115 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3116 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3117 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift);CHKERRQ(ierr); 3118 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3119 diag += 16; 3120 } 3121 break; 3122 case 5: 3123 for (i=0; i<mbs; i++) { 3124 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3125 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3126 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift);CHKERRQ(ierr); 3127 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3128 diag += 25; 3129 } 3130 break; 3131 case 6: 3132 for (i=0; i<mbs; i++) { 3133 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; 3134 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3135 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift);CHKERRQ(ierr); 3136 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3137 diag += 36; 3138 } 3139 break; 3140 case 7: 3141 for (i=0; i<mbs; i++) { 3142 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; 3143 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3144 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift);CHKERRQ(ierr); 3145 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3146 diag += 49; 3147 } 3148 break; 3149 default: 3150 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3151 for (i=0; i<mbs; i++) { 3152 for (j=0; j<bs; j++) { 3153 IJ[j] = bs*i + j; 3154 } 3155 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3156 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);CHKERRQ(ierr); 3157 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3158 diag += bs2; 3159 } 3160 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3161 } 3162 a->ibdiagvalid = PETSC_TRUE; 3163 PetscFunctionReturn(0); 3164 } 3165 3166 #undef __FUNCT__ 3167 #define __FUNCT__ "MatSetRandom_SeqAIJ" 3168 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3169 { 3170 PetscErrorCode ierr; 3171 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3172 PetscScalar a; 3173 PetscInt m,n,i,j,col; 3174 3175 PetscFunctionBegin; 3176 if (!x->assembled) { 3177 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3178 for (i=0; i<m; i++) { 3179 for (j=0; j<aij->imax[i]; j++) { 3180 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3181 col = (PetscInt)(n*PetscRealPart(a)); 3182 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3183 } 3184 } 3185 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3186 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3187 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3188 PetscFunctionReturn(0); 3189 } 3190 3191 /* -------------------------------------------------------------------*/ 3192 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3193 MatGetRow_SeqAIJ, 3194 MatRestoreRow_SeqAIJ, 3195 MatMult_SeqAIJ, 3196 /* 4*/ MatMultAdd_SeqAIJ, 3197 MatMultTranspose_SeqAIJ, 3198 MatMultTransposeAdd_SeqAIJ, 3199 0, 3200 0, 3201 0, 3202 /* 10*/ 0, 3203 MatLUFactor_SeqAIJ, 3204 0, 3205 MatSOR_SeqAIJ, 3206 MatTranspose_SeqAIJ, 3207 /*1 5*/ MatGetInfo_SeqAIJ, 3208 MatEqual_SeqAIJ, 3209 MatGetDiagonal_SeqAIJ, 3210 MatDiagonalScale_SeqAIJ, 3211 MatNorm_SeqAIJ, 3212 /* 20*/ 0, 3213 MatAssemblyEnd_SeqAIJ, 3214 MatSetOption_SeqAIJ, 3215 MatZeroEntries_SeqAIJ, 3216 /* 24*/ MatZeroRows_SeqAIJ, 3217 0, 3218 0, 3219 0, 3220 0, 3221 /* 29*/ MatSetUp_SeqAIJ, 3222 0, 3223 0, 3224 0, 3225 0, 3226 /* 34*/ MatDuplicate_SeqAIJ, 3227 0, 3228 0, 3229 MatILUFactor_SeqAIJ, 3230 0, 3231 /* 39*/ MatAXPY_SeqAIJ, 3232 MatGetSubMatrices_SeqAIJ, 3233 MatIncreaseOverlap_SeqAIJ, 3234 MatGetValues_SeqAIJ, 3235 MatCopy_SeqAIJ, 3236 /* 44*/ MatGetRowMax_SeqAIJ, 3237 MatScale_SeqAIJ, 3238 0, 3239 MatDiagonalSet_SeqAIJ, 3240 MatZeroRowsColumns_SeqAIJ, 3241 /* 49*/ MatSetRandom_SeqAIJ, 3242 MatGetRowIJ_SeqAIJ, 3243 MatRestoreRowIJ_SeqAIJ, 3244 MatGetColumnIJ_SeqAIJ, 3245 MatRestoreColumnIJ_SeqAIJ, 3246 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3247 0, 3248 0, 3249 MatPermute_SeqAIJ, 3250 0, 3251 /* 59*/ 0, 3252 MatDestroy_SeqAIJ, 3253 MatView_SeqAIJ, 3254 0, 3255 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3256 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3257 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3258 0, 3259 0, 3260 0, 3261 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3262 MatGetRowMinAbs_SeqAIJ, 3263 0, 3264 MatSetColoring_SeqAIJ, 3265 0, 3266 /* 74*/ MatSetValuesAdifor_SeqAIJ, 3267 MatFDColoringApply_AIJ, 3268 0, 3269 0, 3270 0, 3271 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3272 0, 3273 0, 3274 0, 3275 MatLoad_SeqAIJ, 3276 /* 84*/ MatIsSymmetric_SeqAIJ, 3277 MatIsHermitian_SeqAIJ, 3278 0, 3279 0, 3280 0, 3281 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3282 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3283 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3284 MatPtAP_SeqAIJ_SeqAIJ, 3285 MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy, 3286 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3287 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3288 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3289 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3290 0, 3291 /* 99*/ 0, 3292 0, 3293 0, 3294 MatConjugate_SeqAIJ, 3295 0, 3296 /*104*/ MatSetValuesRow_SeqAIJ, 3297 MatRealPart_SeqAIJ, 3298 MatImaginaryPart_SeqAIJ, 3299 0, 3300 0, 3301 /*109*/ MatMatSolve_SeqAIJ, 3302 0, 3303 MatGetRowMin_SeqAIJ, 3304 0, 3305 MatMissingDiagonal_SeqAIJ, 3306 /*114*/ 0, 3307 0, 3308 0, 3309 0, 3310 0, 3311 /*119*/ 0, 3312 0, 3313 0, 3314 0, 3315 MatGetMultiProcBlock_SeqAIJ, 3316 /*124*/ MatFindNonzeroRows_SeqAIJ, 3317 MatGetColumnNorms_SeqAIJ, 3318 MatInvertBlockDiagonal_SeqAIJ, 3319 0, 3320 0, 3321 /*129*/ 0, 3322 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3323 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3324 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3325 MatTransposeColoringCreate_SeqAIJ, 3326 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3327 MatTransColoringApplyDenToSp_SeqAIJ, 3328 MatRARt_SeqAIJ_SeqAIJ, 3329 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3330 MatRARtNumeric_SeqAIJ_SeqAIJ, 3331 /*139*/0, 3332 0, 3333 0, 3334 MatFDColoringSetUp_SeqXAIJ, 3335 MatFindOffBlockDiagonalEntries_SeqAIJ, 3336 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ 3337 }; 3338 3339 #undef __FUNCT__ 3340 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3341 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3342 { 3343 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3344 PetscInt i,nz,n; 3345 3346 PetscFunctionBegin; 3347 nz = aij->maxnz; 3348 n = mat->rmap->n; 3349 for (i=0; i<nz; i++) { 3350 aij->j[i] = indices[i]; 3351 } 3352 aij->nz = nz; 3353 for (i=0; i<n; i++) { 3354 aij->ilen[i] = aij->imax[i]; 3355 } 3356 PetscFunctionReturn(0); 3357 } 3358 3359 #undef __FUNCT__ 3360 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3361 /*@ 3362 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3363 in the matrix. 3364 3365 Input Parameters: 3366 + mat - the SeqAIJ matrix 3367 - indices - the column indices 3368 3369 Level: advanced 3370 3371 Notes: 3372 This can be called if you have precomputed the nonzero structure of the 3373 matrix and want to provide it to the matrix object to improve the performance 3374 of the MatSetValues() operation. 3375 3376 You MUST have set the correct numbers of nonzeros per row in the call to 3377 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3378 3379 MUST be called before any calls to MatSetValues(); 3380 3381 The indices should start with zero, not one. 3382 3383 @*/ 3384 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3385 { 3386 PetscErrorCode ierr; 3387 3388 PetscFunctionBegin; 3389 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3390 PetscValidPointer(indices,2); 3391 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3392 PetscFunctionReturn(0); 3393 } 3394 3395 /* ----------------------------------------------------------------------------------------*/ 3396 3397 #undef __FUNCT__ 3398 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3399 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3400 { 3401 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3402 PetscErrorCode ierr; 3403 size_t nz = aij->i[mat->rmap->n]; 3404 3405 PetscFunctionBegin; 3406 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3407 3408 /* allocate space for values if not already there */ 3409 if (!aij->saved_values) { 3410 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3411 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3412 } 3413 3414 /* copy values over */ 3415 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3416 PetscFunctionReturn(0); 3417 } 3418 3419 #undef __FUNCT__ 3420 #define __FUNCT__ "MatStoreValues" 3421 /*@ 3422 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3423 example, reuse of the linear part of a Jacobian, while recomputing the 3424 nonlinear portion. 3425 3426 Collect on Mat 3427 3428 Input Parameters: 3429 . mat - the matrix (currently only AIJ matrices support this option) 3430 3431 Level: advanced 3432 3433 Common Usage, with SNESSolve(): 3434 $ Create Jacobian matrix 3435 $ Set linear terms into matrix 3436 $ Apply boundary conditions to matrix, at this time matrix must have 3437 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3438 $ boundary conditions again will not change the nonzero structure 3439 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3440 $ ierr = MatStoreValues(mat); 3441 $ Call SNESSetJacobian() with matrix 3442 $ In your Jacobian routine 3443 $ ierr = MatRetrieveValues(mat); 3444 $ Set nonlinear terms in matrix 3445 3446 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3447 $ // build linear portion of Jacobian 3448 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3449 $ ierr = MatStoreValues(mat); 3450 $ loop over nonlinear iterations 3451 $ ierr = MatRetrieveValues(mat); 3452 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3453 $ // call MatAssemblyBegin/End() on matrix 3454 $ Solve linear system with Jacobian 3455 $ endloop 3456 3457 Notes: 3458 Matrix must already be assemblied before calling this routine 3459 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3460 calling this routine. 3461 3462 When this is called multiple times it overwrites the previous set of stored values 3463 and does not allocated additional space. 3464 3465 .seealso: MatRetrieveValues() 3466 3467 @*/ 3468 PetscErrorCode MatStoreValues(Mat mat) 3469 { 3470 PetscErrorCode ierr; 3471 3472 PetscFunctionBegin; 3473 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3474 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3475 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3476 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3477 PetscFunctionReturn(0); 3478 } 3479 3480 #undef __FUNCT__ 3481 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3482 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3483 { 3484 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3485 PetscErrorCode ierr; 3486 PetscInt nz = aij->i[mat->rmap->n]; 3487 3488 PetscFunctionBegin; 3489 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3490 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3491 /* copy values over */ 3492 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3493 PetscFunctionReturn(0); 3494 } 3495 3496 #undef __FUNCT__ 3497 #define __FUNCT__ "MatRetrieveValues" 3498 /*@ 3499 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3500 example, reuse of the linear part of a Jacobian, while recomputing the 3501 nonlinear portion. 3502 3503 Collect on Mat 3504 3505 Input Parameters: 3506 . mat - the matrix (currently on AIJ matrices support this option) 3507 3508 Level: advanced 3509 3510 .seealso: MatStoreValues() 3511 3512 @*/ 3513 PetscErrorCode MatRetrieveValues(Mat mat) 3514 { 3515 PetscErrorCode ierr; 3516 3517 PetscFunctionBegin; 3518 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3519 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3520 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3521 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3522 PetscFunctionReturn(0); 3523 } 3524 3525 3526 /* --------------------------------------------------------------------------------*/ 3527 #undef __FUNCT__ 3528 #define __FUNCT__ "MatCreateSeqAIJ" 3529 /*@C 3530 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3531 (the default parallel PETSc format). For good matrix assembly performance 3532 the user should preallocate the matrix storage by setting the parameter nz 3533 (or the array nnz). By setting these parameters accurately, performance 3534 during matrix assembly can be increased by more than a factor of 50. 3535 3536 Collective on MPI_Comm 3537 3538 Input Parameters: 3539 + comm - MPI communicator, set to PETSC_COMM_SELF 3540 . m - number of rows 3541 . n - number of columns 3542 . nz - number of nonzeros per row (same for all rows) 3543 - nnz - array containing the number of nonzeros in the various rows 3544 (possibly different for each row) or NULL 3545 3546 Output Parameter: 3547 . A - the matrix 3548 3549 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3550 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3551 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3552 3553 Notes: 3554 If nnz is given then nz is ignored 3555 3556 The AIJ format (also called the Yale sparse matrix format or 3557 compressed row storage), is fully compatible with standard Fortran 77 3558 storage. That is, the stored row and column indices can begin at 3559 either one (as in Fortran) or zero. See the users' manual for details. 3560 3561 Specify the preallocated storage with either nz or nnz (not both). 3562 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3563 allocation. For large problems you MUST preallocate memory or you 3564 will get TERRIBLE performance, see the users' manual chapter on matrices. 3565 3566 By default, this format uses inodes (identical nodes) when possible, to 3567 improve numerical efficiency of matrix-vector products and solves. We 3568 search for consecutive rows with the same nonzero structure, thereby 3569 reusing matrix information to achieve increased efficiency. 3570 3571 Options Database Keys: 3572 + -mat_no_inode - Do not use inodes 3573 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3574 3575 Level: intermediate 3576 3577 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3578 3579 @*/ 3580 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3581 { 3582 PetscErrorCode ierr; 3583 3584 PetscFunctionBegin; 3585 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3586 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3587 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3588 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3589 PetscFunctionReturn(0); 3590 } 3591 3592 #undef __FUNCT__ 3593 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3594 /*@C 3595 MatSeqAIJSetPreallocation - For good matrix assembly performance 3596 the user should preallocate the matrix storage by setting the parameter nz 3597 (or the array nnz). By setting these parameters accurately, performance 3598 during matrix assembly can be increased by more than a factor of 50. 3599 3600 Collective on MPI_Comm 3601 3602 Input Parameters: 3603 + B - The matrix 3604 . nz - number of nonzeros per row (same for all rows) 3605 - nnz - array containing the number of nonzeros in the various rows 3606 (possibly different for each row) or NULL 3607 3608 Notes: 3609 If nnz is given then nz is ignored 3610 3611 The AIJ format (also called the Yale sparse matrix format or 3612 compressed row storage), is fully compatible with standard Fortran 77 3613 storage. That is, the stored row and column indices can begin at 3614 either one (as in Fortran) or zero. See the users' manual for details. 3615 3616 Specify the preallocated storage with either nz or nnz (not both). 3617 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3618 allocation. For large problems you MUST preallocate memory or you 3619 will get TERRIBLE performance, see the users' manual chapter on matrices. 3620 3621 You can call MatGetInfo() to get information on how effective the preallocation was; 3622 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3623 You can also run with the option -info and look for messages with the string 3624 malloc in them to see if additional memory allocation was needed. 3625 3626 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3627 entries or columns indices 3628 3629 By default, this format uses inodes (identical nodes) when possible, to 3630 improve numerical efficiency of matrix-vector products and solves. We 3631 search for consecutive rows with the same nonzero structure, thereby 3632 reusing matrix information to achieve increased efficiency. 3633 3634 Options Database Keys: 3635 + -mat_no_inode - Do not use inodes 3636 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3637 - -mat_aij_oneindex - Internally use indexing starting at 1 3638 rather than 0. Note that when calling MatSetValues(), 3639 the user still MUST index entries starting at 0! 3640 3641 Level: intermediate 3642 3643 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3644 3645 @*/ 3646 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3647 { 3648 PetscErrorCode ierr; 3649 3650 PetscFunctionBegin; 3651 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3652 PetscValidType(B,1); 3653 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3654 PetscFunctionReturn(0); 3655 } 3656 3657 #undef __FUNCT__ 3658 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3659 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3660 { 3661 Mat_SeqAIJ *b; 3662 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3663 PetscErrorCode ierr; 3664 PetscInt i; 3665 3666 PetscFunctionBegin; 3667 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3668 if (nz == MAT_SKIP_ALLOCATION) { 3669 skipallocation = PETSC_TRUE; 3670 nz = 0; 3671 } 3672 3673 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3674 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3675 3676 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3677 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3678 if (nnz) { 3679 for (i=0; i<B->rmap->n; i++) { 3680 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]); 3681 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); 3682 } 3683 } 3684 3685 B->preallocated = PETSC_TRUE; 3686 3687 b = (Mat_SeqAIJ*)B->data; 3688 3689 if (!skipallocation) { 3690 if (!b->imax) { 3691 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3692 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3693 } 3694 if (!nnz) { 3695 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3696 else if (nz < 0) nz = 1; 3697 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3698 nz = nz*B->rmap->n; 3699 } else { 3700 nz = 0; 3701 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3702 } 3703 /* b->ilen will count nonzeros in each row so far. */ 3704 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3705 3706 /* allocate the matrix space */ 3707 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3708 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3709 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3710 b->i[0] = 0; 3711 for (i=1; i<B->rmap->n+1; i++) { 3712 b->i[i] = b->i[i-1] + b->imax[i-1]; 3713 } 3714 b->singlemalloc = PETSC_TRUE; 3715 b->free_a = PETSC_TRUE; 3716 b->free_ij = PETSC_TRUE; 3717 #if defined(PETSC_THREADCOMM_ACTIVE) 3718 ierr = MatZeroEntries_SeqAIJ(B);CHKERRQ(ierr); 3719 #endif 3720 } else { 3721 b->free_a = PETSC_FALSE; 3722 b->free_ij = PETSC_FALSE; 3723 } 3724 3725 b->nz = 0; 3726 b->maxnz = nz; 3727 B->info.nz_unneeded = (double)b->maxnz; 3728 if (realalloc) { 3729 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3730 } 3731 PetscFunctionReturn(0); 3732 } 3733 3734 #undef __FUNCT__ 3735 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3736 /*@ 3737 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3738 3739 Input Parameters: 3740 + B - the matrix 3741 . i - the indices into j for the start of each row (starts with zero) 3742 . j - the column indices for each row (starts with zero) these must be sorted for each row 3743 - v - optional values in the matrix 3744 3745 Level: developer 3746 3747 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3748 3749 .keywords: matrix, aij, compressed row, sparse, sequential 3750 3751 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3752 @*/ 3753 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3754 { 3755 PetscErrorCode ierr; 3756 3757 PetscFunctionBegin; 3758 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3759 PetscValidType(B,1); 3760 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3761 PetscFunctionReturn(0); 3762 } 3763 3764 #undef __FUNCT__ 3765 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3766 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3767 { 3768 PetscInt i; 3769 PetscInt m,n; 3770 PetscInt nz; 3771 PetscInt *nnz, nz_max = 0; 3772 PetscScalar *values; 3773 PetscErrorCode ierr; 3774 3775 PetscFunctionBegin; 3776 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3777 3778 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3779 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3780 3781 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3782 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3783 for (i = 0; i < m; i++) { 3784 nz = Ii[i+1]- Ii[i]; 3785 nz_max = PetscMax(nz_max, nz); 3786 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3787 nnz[i] = nz; 3788 } 3789 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3790 ierr = PetscFree(nnz);CHKERRQ(ierr); 3791 3792 if (v) { 3793 values = (PetscScalar*) v; 3794 } else { 3795 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3796 } 3797 3798 for (i = 0; i < m; i++) { 3799 nz = Ii[i+1] - Ii[i]; 3800 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3801 } 3802 3803 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3804 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3805 3806 if (!v) { 3807 ierr = PetscFree(values);CHKERRQ(ierr); 3808 } 3809 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3810 PetscFunctionReturn(0); 3811 } 3812 3813 #include <../src/mat/impls/dense/seq/dense.h> 3814 #include <petsc-private/kernels/petscaxpy.h> 3815 3816 #undef __FUNCT__ 3817 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3818 /* 3819 Computes (B'*A')' since computing B*A directly is untenable 3820 3821 n p p 3822 ( ) ( ) ( ) 3823 m ( A ) * n ( B ) = m ( C ) 3824 ( ) ( ) ( ) 3825 3826 */ 3827 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3828 { 3829 PetscErrorCode ierr; 3830 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3831 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3832 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3833 PetscInt i,n,m,q,p; 3834 const PetscInt *ii,*idx; 3835 const PetscScalar *b,*a,*a_q; 3836 PetscScalar *c,*c_q; 3837 3838 PetscFunctionBegin; 3839 m = A->rmap->n; 3840 n = A->cmap->n; 3841 p = B->cmap->n; 3842 a = sub_a->v; 3843 b = sub_b->a; 3844 c = sub_c->v; 3845 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3846 3847 ii = sub_b->i; 3848 idx = sub_b->j; 3849 for (i=0; i<n; i++) { 3850 q = ii[i+1] - ii[i]; 3851 while (q-->0) { 3852 c_q = c + m*(*idx); 3853 a_q = a + m*i; 3854 PetscKernelAXPY(c_q,*b,a_q,m); 3855 idx++; 3856 b++; 3857 } 3858 } 3859 PetscFunctionReturn(0); 3860 } 3861 3862 #undef __FUNCT__ 3863 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3864 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3865 { 3866 PetscErrorCode ierr; 3867 PetscInt m=A->rmap->n,n=B->cmap->n; 3868 Mat Cmat; 3869 3870 PetscFunctionBegin; 3871 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); 3872 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3873 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3874 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3875 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3876 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3877 3878 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3879 3880 *C = Cmat; 3881 PetscFunctionReturn(0); 3882 } 3883 3884 /* ----------------------------------------------------------------*/ 3885 #undef __FUNCT__ 3886 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3887 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3888 { 3889 PetscErrorCode ierr; 3890 3891 PetscFunctionBegin; 3892 if (scall == MAT_INITIAL_MATRIX) { 3893 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3894 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3895 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3896 } 3897 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3898 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3899 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3900 PetscFunctionReturn(0); 3901 } 3902 3903 3904 /*MC 3905 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3906 based on compressed sparse row format. 3907 3908 Options Database Keys: 3909 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3910 3911 Level: beginner 3912 3913 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3914 M*/ 3915 3916 /*MC 3917 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3918 3919 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3920 and MATMPIAIJ otherwise. As a result, for single process communicators, 3921 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3922 for communicators controlling multiple processes. It is recommended that you call both of 3923 the above preallocation routines for simplicity. 3924 3925 Options Database Keys: 3926 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3927 3928 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3929 enough exist. 3930 3931 Level: beginner 3932 3933 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3934 M*/ 3935 3936 /*MC 3937 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3938 3939 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3940 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3941 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3942 for communicators controlling multiple processes. It is recommended that you call both of 3943 the above preallocation routines for simplicity. 3944 3945 Options Database Keys: 3946 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3947 3948 Level: beginner 3949 3950 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3951 M*/ 3952 3953 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3954 #if defined(PETSC_HAVE_ELEMENTAL) 3955 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3956 #endif 3957 3958 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3959 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3960 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3961 #endif 3962 3963 3964 #undef __FUNCT__ 3965 #define __FUNCT__ "MatSeqAIJGetArray" 3966 /*@C 3967 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 3968 3969 Not Collective 3970 3971 Input Parameter: 3972 . mat - a MATSEQAIJ matrix 3973 3974 Output Parameter: 3975 . array - pointer to the data 3976 3977 Level: intermediate 3978 3979 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3980 @*/ 3981 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3982 { 3983 PetscErrorCode ierr; 3984 3985 PetscFunctionBegin; 3986 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3987 PetscFunctionReturn(0); 3988 } 3989 3990 #undef __FUNCT__ 3991 #define __FUNCT__ "MatSeqAIJGetMaxRowNonzeros" 3992 /*@C 3993 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 3994 3995 Not Collective 3996 3997 Input Parameter: 3998 . mat - a MATSEQAIJ matrix 3999 4000 Output Parameter: 4001 . nz - the maximum number of nonzeros in any row 4002 4003 Level: intermediate 4004 4005 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4006 @*/ 4007 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4008 { 4009 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4010 4011 PetscFunctionBegin; 4012 *nz = aij->rmax; 4013 PetscFunctionReturn(0); 4014 } 4015 4016 #undef __FUNCT__ 4017 #define __FUNCT__ "MatSeqAIJRestoreArray" 4018 /*@C 4019 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4020 4021 Not Collective 4022 4023 Input Parameters: 4024 . mat - a MATSEQAIJ matrix 4025 . array - pointer to the data 4026 4027 Level: intermediate 4028 4029 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4030 @*/ 4031 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4032 { 4033 PetscErrorCode ierr; 4034 4035 PetscFunctionBegin; 4036 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4037 PetscFunctionReturn(0); 4038 } 4039 4040 #undef __FUNCT__ 4041 #define __FUNCT__ "MatCreate_SeqAIJ" 4042 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4043 { 4044 Mat_SeqAIJ *b; 4045 PetscErrorCode ierr; 4046 PetscMPIInt size; 4047 4048 PetscFunctionBegin; 4049 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4050 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4051 4052 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4053 4054 B->data = (void*)b; 4055 4056 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4057 4058 b->row = 0; 4059 b->col = 0; 4060 b->icol = 0; 4061 b->reallocs = 0; 4062 b->ignorezeroentries = PETSC_FALSE; 4063 b->roworiented = PETSC_TRUE; 4064 b->nonew = 0; 4065 b->diag = 0; 4066 b->solve_work = 0; 4067 B->spptr = 0; 4068 b->saved_values = 0; 4069 b->idiag = 0; 4070 b->mdiag = 0; 4071 b->ssor_work = 0; 4072 b->omega = 1.0; 4073 b->fshift = 0.0; 4074 b->idiagvalid = PETSC_FALSE; 4075 b->ibdiagvalid = PETSC_FALSE; 4076 b->keepnonzeropattern = PETSC_FALSE; 4077 4078 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4079 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4080 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4081 4082 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4083 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4084 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4085 #endif 4086 4087 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4088 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4089 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4090 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4091 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4092 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4093 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4094 #if defined(PETSC_HAVE_ELEMENTAL) 4095 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4096 #endif 4097 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4098 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4099 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4100 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4101 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4102 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4103 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4104 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4105 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4106 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4107 PetscFunctionReturn(0); 4108 } 4109 4110 #undef __FUNCT__ 4111 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 4112 /* 4113 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4114 */ 4115 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4116 { 4117 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4118 PetscErrorCode ierr; 4119 PetscInt i,m = A->rmap->n; 4120 4121 PetscFunctionBegin; 4122 c = (Mat_SeqAIJ*)C->data; 4123 4124 C->factortype = A->factortype; 4125 c->row = 0; 4126 c->col = 0; 4127 c->icol = 0; 4128 c->reallocs = 0; 4129 4130 C->assembled = PETSC_TRUE; 4131 4132 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4133 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4134 4135 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4136 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4137 for (i=0; i<m; i++) { 4138 c->imax[i] = a->imax[i]; 4139 c->ilen[i] = a->ilen[i]; 4140 } 4141 4142 /* allocate the matrix space */ 4143 if (mallocmatspace) { 4144 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4145 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4146 4147 c->singlemalloc = PETSC_TRUE; 4148 4149 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4150 if (m > 0) { 4151 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4152 if (cpvalues == MAT_COPY_VALUES) { 4153 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4154 } else { 4155 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4156 } 4157 } 4158 } 4159 4160 c->ignorezeroentries = a->ignorezeroentries; 4161 c->roworiented = a->roworiented; 4162 c->nonew = a->nonew; 4163 if (a->diag) { 4164 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4165 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4166 for (i=0; i<m; i++) { 4167 c->diag[i] = a->diag[i]; 4168 } 4169 } else c->diag = 0; 4170 4171 c->solve_work = 0; 4172 c->saved_values = 0; 4173 c->idiag = 0; 4174 c->ssor_work = 0; 4175 c->keepnonzeropattern = a->keepnonzeropattern; 4176 c->free_a = PETSC_TRUE; 4177 c->free_ij = PETSC_TRUE; 4178 4179 c->rmax = a->rmax; 4180 c->nz = a->nz; 4181 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4182 C->preallocated = PETSC_TRUE; 4183 4184 c->compressedrow.use = a->compressedrow.use; 4185 c->compressedrow.nrows = a->compressedrow.nrows; 4186 if (a->compressedrow.use) { 4187 i = a->compressedrow.nrows; 4188 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4189 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4190 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4191 } else { 4192 c->compressedrow.use = PETSC_FALSE; 4193 c->compressedrow.i = NULL; 4194 c->compressedrow.rindex = NULL; 4195 } 4196 C->nonzerostate = A->nonzerostate; 4197 4198 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4199 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4200 PetscFunctionReturn(0); 4201 } 4202 4203 #undef __FUNCT__ 4204 #define __FUNCT__ "MatDuplicate_SeqAIJ" 4205 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4206 { 4207 PetscErrorCode ierr; 4208 4209 PetscFunctionBegin; 4210 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4211 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4212 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4213 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4214 } 4215 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4216 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4217 PetscFunctionReturn(0); 4218 } 4219 4220 #undef __FUNCT__ 4221 #define __FUNCT__ "MatLoad_SeqAIJ" 4222 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4223 { 4224 Mat_SeqAIJ *a; 4225 PetscErrorCode ierr; 4226 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4227 int fd; 4228 PetscMPIInt size; 4229 MPI_Comm comm; 4230 PetscInt bs = newMat->rmap->bs; 4231 4232 PetscFunctionBegin; 4233 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4234 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4235 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4236 4237 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4238 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4239 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4240 if (bs < 0) bs = 1; 4241 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4242 4243 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4244 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4245 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4246 M = header[1]; N = header[2]; nz = header[3]; 4247 4248 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4249 4250 /* read in row lengths */ 4251 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4252 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4253 4254 /* check if sum of rowlengths is same as nz */ 4255 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4256 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); 4257 4258 /* set global size if not set already*/ 4259 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4260 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4261 } else { 4262 /* if sizes and type are already set, check if the vector global sizes are correct */ 4263 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4264 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4265 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4266 } 4267 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); 4268 } 4269 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4270 a = (Mat_SeqAIJ*)newMat->data; 4271 4272 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4273 4274 /* read in nonzero values */ 4275 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4276 4277 /* set matrix "i" values */ 4278 a->i[0] = 0; 4279 for (i=1; i<= M; i++) { 4280 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4281 a->ilen[i-1] = rowlengths[i-1]; 4282 } 4283 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4284 4285 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4286 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4287 PetscFunctionReturn(0); 4288 } 4289 4290 #undef __FUNCT__ 4291 #define __FUNCT__ "MatEqual_SeqAIJ" 4292 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4293 { 4294 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4295 PetscErrorCode ierr; 4296 #if defined(PETSC_USE_COMPLEX) 4297 PetscInt k; 4298 #endif 4299 4300 PetscFunctionBegin; 4301 /* If the matrix dimensions are not equal,or no of nonzeros */ 4302 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4303 *flg = PETSC_FALSE; 4304 PetscFunctionReturn(0); 4305 } 4306 4307 /* if the a->i are the same */ 4308 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4309 if (!*flg) PetscFunctionReturn(0); 4310 4311 /* if a->j are the same */ 4312 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4313 if (!*flg) PetscFunctionReturn(0); 4314 4315 /* if a->a are the same */ 4316 #if defined(PETSC_USE_COMPLEX) 4317 for (k=0; k<a->nz; k++) { 4318 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4319 *flg = PETSC_FALSE; 4320 PetscFunctionReturn(0); 4321 } 4322 } 4323 #else 4324 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4325 #endif 4326 PetscFunctionReturn(0); 4327 } 4328 4329 #undef __FUNCT__ 4330 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4331 /*@ 4332 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4333 provided by the user. 4334 4335 Collective on MPI_Comm 4336 4337 Input Parameters: 4338 + comm - must be an MPI communicator of size 1 4339 . m - number of rows 4340 . n - number of columns 4341 . i - row indices 4342 . j - column indices 4343 - a - matrix values 4344 4345 Output Parameter: 4346 . mat - the matrix 4347 4348 Level: intermediate 4349 4350 Notes: 4351 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4352 once the matrix is destroyed and not before 4353 4354 You cannot set new nonzero locations into this matrix, that will generate an error. 4355 4356 The i and j indices are 0 based 4357 4358 The format which is used for the sparse matrix input, is equivalent to a 4359 row-major ordering.. i.e for the following matrix, the input data expected is 4360 as shown: 4361 4362 1 0 0 4363 2 0 3 4364 4 5 6 4365 4366 i = {0,1,3,6} [size = nrow+1 = 3+1] 4367 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4368 v = {1,2,3,4,5,6} [size = nz = 6] 4369 4370 4371 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4372 4373 @*/ 4374 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 4375 { 4376 PetscErrorCode ierr; 4377 PetscInt ii; 4378 Mat_SeqAIJ *aij; 4379 #if defined(PETSC_USE_DEBUG) 4380 PetscInt jj; 4381 #endif 4382 4383 PetscFunctionBegin; 4384 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4385 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4386 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4387 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4388 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4389 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4390 aij = (Mat_SeqAIJ*)(*mat)->data; 4391 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4392 4393 aij->i = i; 4394 aij->j = j; 4395 aij->a = a; 4396 aij->singlemalloc = PETSC_FALSE; 4397 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4398 aij->free_a = PETSC_FALSE; 4399 aij->free_ij = PETSC_FALSE; 4400 4401 for (ii=0; ii<m; ii++) { 4402 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4403 #if defined(PETSC_USE_DEBUG) 4404 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]); 4405 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4406 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); 4407 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); 4408 } 4409 #endif 4410 } 4411 #if defined(PETSC_USE_DEBUG) 4412 for (ii=0; ii<aij->i[m]; ii++) { 4413 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4414 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]); 4415 } 4416 #endif 4417 4418 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4419 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4420 PetscFunctionReturn(0); 4421 } 4422 #undef __FUNCT__ 4423 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4424 /*@C 4425 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4426 provided by the user. 4427 4428 Collective on MPI_Comm 4429 4430 Input Parameters: 4431 + comm - must be an MPI communicator of size 1 4432 . m - number of rows 4433 . n - number of columns 4434 . i - row indices 4435 . j - column indices 4436 . a - matrix values 4437 . nz - number of nonzeros 4438 - idx - 0 or 1 based 4439 4440 Output Parameter: 4441 . mat - the matrix 4442 4443 Level: intermediate 4444 4445 Notes: 4446 The i and j indices are 0 based 4447 4448 The format which is used for the sparse matrix input, is equivalent to a 4449 row-major ordering.. i.e for the following matrix, the input data expected is 4450 as shown: 4451 4452 1 0 0 4453 2 0 3 4454 4 5 6 4455 4456 i = {0,1,1,2,2,2} 4457 j = {0,0,2,0,1,2} 4458 v = {1,2,3,4,5,6} 4459 4460 4461 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4462 4463 @*/ 4464 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4465 { 4466 PetscErrorCode ierr; 4467 PetscInt ii, *nnz, one = 1,row,col; 4468 4469 4470 PetscFunctionBegin; 4471 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4472 for (ii = 0; ii < nz; ii++) { 4473 nnz[i[ii] - !!idx] += 1; 4474 } 4475 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4476 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4477 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4478 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4479 for (ii = 0; ii < nz; ii++) { 4480 if (idx) { 4481 row = i[ii] - 1; 4482 col = j[ii] - 1; 4483 } else { 4484 row = i[ii]; 4485 col = j[ii]; 4486 } 4487 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4488 } 4489 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4490 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4491 ierr = PetscFree(nnz);CHKERRQ(ierr); 4492 PetscFunctionReturn(0); 4493 } 4494 4495 #undef __FUNCT__ 4496 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4497 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4498 { 4499 PetscErrorCode ierr; 4500 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4501 4502 PetscFunctionBegin; 4503 if (coloring->ctype == IS_COLORING_GLOBAL) { 4504 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4505 a->coloring = coloring; 4506 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4507 PetscInt i,*larray; 4508 ISColoring ocoloring; 4509 ISColoringValue *colors; 4510 4511 /* set coloring for diagonal portion */ 4512 ierr = PetscMalloc1(A->cmap->n,&larray);CHKERRQ(ierr); 4513 for (i=0; i<A->cmap->n; i++) larray[i] = i; 4514 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,NULL,larray);CHKERRQ(ierr); 4515 ierr = PetscMalloc1(A->cmap->n,&colors);CHKERRQ(ierr); 4516 for (i=0; i<A->cmap->n; i++) colors[i] = coloring->colors[larray[i]]; 4517 ierr = PetscFree(larray);CHKERRQ(ierr); 4518 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4519 a->coloring = ocoloring; 4520 } 4521 PetscFunctionReturn(0); 4522 } 4523 4524 #undef __FUNCT__ 4525 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4526 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4527 { 4528 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4529 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4530 MatScalar *v = a->a; 4531 PetscScalar *values = (PetscScalar*)advalues; 4532 ISColoringValue *color; 4533 4534 PetscFunctionBegin; 4535 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4536 color = a->coloring->colors; 4537 /* loop over rows */ 4538 for (i=0; i<m; i++) { 4539 nz = ii[i+1] - ii[i]; 4540 /* loop over columns putting computed value into matrix */ 4541 for (j=0; j<nz; j++) *v++ = values[color[*jj++]]; 4542 values += nl; /* jump to next row of derivatives */ 4543 } 4544 PetscFunctionReturn(0); 4545 } 4546 4547 #undef __FUNCT__ 4548 #define __FUNCT__ "MatSeqAIJInvalidateDiagonal" 4549 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4550 { 4551 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4552 PetscErrorCode ierr; 4553 4554 PetscFunctionBegin; 4555 a->idiagvalid = PETSC_FALSE; 4556 a->ibdiagvalid = PETSC_FALSE; 4557 4558 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4559 PetscFunctionReturn(0); 4560 } 4561 4562 #undef __FUNCT__ 4563 #define __FUNCT__ "MatCreateMPIMatConcatenateSeqMat_SeqAIJ" 4564 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4565 { 4566 PetscErrorCode ierr; 4567 4568 PetscFunctionBegin; 4569 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4570 PetscFunctionReturn(0); 4571 } 4572 4573 /* 4574 Special version for direct calls from Fortran 4575 */ 4576 #include <petsc-private/fortranimpl.h> 4577 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4578 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4579 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4580 #define matsetvaluesseqaij_ matsetvaluesseqaij 4581 #endif 4582 4583 /* Change these macros so can be used in void function */ 4584 #undef CHKERRQ 4585 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4586 #undef SETERRQ2 4587 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4588 #undef SETERRQ3 4589 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4590 4591 #undef __FUNCT__ 4592 #define __FUNCT__ "matsetvaluesseqaij_" 4593 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) 4594 { 4595 Mat A = *AA; 4596 PetscInt m = *mm, n = *nn; 4597 InsertMode is = *isis; 4598 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4599 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4600 PetscInt *imax,*ai,*ailen; 4601 PetscErrorCode ierr; 4602 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4603 MatScalar *ap,value,*aa; 4604 PetscBool ignorezeroentries = a->ignorezeroentries; 4605 PetscBool roworiented = a->roworiented; 4606 4607 PetscFunctionBegin; 4608 MatCheckPreallocated(A,1); 4609 imax = a->imax; 4610 ai = a->i; 4611 ailen = a->ilen; 4612 aj = a->j; 4613 aa = a->a; 4614 4615 for (k=0; k<m; k++) { /* loop over added rows */ 4616 row = im[k]; 4617 if (row < 0) continue; 4618 #if defined(PETSC_USE_DEBUG) 4619 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4620 #endif 4621 rp = aj + ai[row]; ap = aa + ai[row]; 4622 rmax = imax[row]; nrow = ailen[row]; 4623 low = 0; 4624 high = nrow; 4625 for (l=0; l<n; l++) { /* loop over added columns */ 4626 if (in[l] < 0) continue; 4627 #if defined(PETSC_USE_DEBUG) 4628 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4629 #endif 4630 col = in[l]; 4631 if (roworiented) value = v[l + k*n]; 4632 else value = v[k + l*m]; 4633 4634 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4635 4636 if (col <= lastcol) low = 0; 4637 else high = nrow; 4638 lastcol = col; 4639 while (high-low > 5) { 4640 t = (low+high)/2; 4641 if (rp[t] > col) high = t; 4642 else low = t; 4643 } 4644 for (i=low; i<high; i++) { 4645 if (rp[i] > col) break; 4646 if (rp[i] == col) { 4647 if (is == ADD_VALUES) ap[i] += value; 4648 else ap[i] = value; 4649 goto noinsert; 4650 } 4651 } 4652 if (value == 0.0 && ignorezeroentries) goto noinsert; 4653 if (nonew == 1) goto noinsert; 4654 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4655 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4656 N = nrow++ - 1; a->nz++; high++; 4657 /* shift up all the later entries in this row */ 4658 for (ii=N; ii>=i; ii--) { 4659 rp[ii+1] = rp[ii]; 4660 ap[ii+1] = ap[ii]; 4661 } 4662 rp[i] = col; 4663 ap[i] = value; 4664 A->nonzerostate++; 4665 noinsert:; 4666 low = i + 1; 4667 } 4668 ailen[row] = nrow; 4669 } 4670 PetscFunctionReturnVoid(); 4671 } 4672 4673