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