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