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