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