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