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__ "MatGetColumnNorms_SeqAIJ" 15 PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms) 16 { 17 PetscErrorCode ierr; 18 PetscInt i,m,n; 19 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 20 21 PetscFunctionBegin; 22 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 23 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 24 if (type == NORM_2) { 25 for (i=0; i<aij->i[m]; i++) { 26 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 27 } 28 } else if (type == NORM_1) { 29 for (i=0; i<aij->i[m]; i++) { 30 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]); 31 } 32 } else if (type == NORM_INFINITY) { 33 for (i=0; i<aij->i[m]; i++) { 34 norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]); 35 } 36 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 37 38 if (type == NORM_2) { 39 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 40 } 41 PetscFunctionReturn(0); 42 } 43 44 #undef __FUNCT__ 45 #define __FUNCT__ "MatFindOffBlockDiagonalEntries_SeqAIJ" 46 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is) 47 { 48 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 49 PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs; 50 const PetscInt *jj = a->j,*ii = a->i; 51 PetscInt *rows; 52 PetscErrorCode ierr; 53 54 PetscFunctionBegin; 55 for (i=0; i<m; i++) { 56 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 57 cnt++; 58 } 59 } 60 ierr = PetscMalloc1(cnt,&rows);CHKERRQ(ierr); 61 cnt = 0; 62 for (i=0; i<m; i++) { 63 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 64 rows[cnt] = i; 65 cnt++; 66 } 67 } 68 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "MatFindZeroDiagonals_SeqAIJ_Private" 74 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 75 { 76 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 77 const MatScalar *aa = a->a; 78 PetscInt i,m=A->rmap->n,cnt = 0; 79 const PetscInt *ii = a->i,*jj = a->j,*diag; 80 PetscInt *rows; 81 PetscErrorCode ierr; 82 83 PetscFunctionBegin; 84 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 85 diag = a->diag; 86 for (i=0; i<m; i++) { 87 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 88 cnt++; 89 } 90 } 91 ierr = PetscMalloc1(cnt,&rows);CHKERRQ(ierr); 92 cnt = 0; 93 for (i=0; i<m; i++) { 94 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 95 rows[cnt++] = i; 96 } 97 } 98 *nrows = cnt; 99 *zrows = rows; 100 PetscFunctionReturn(0); 101 } 102 103 #undef __FUNCT__ 104 #define __FUNCT__ "MatFindZeroDiagonals_SeqAIJ" 105 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 106 { 107 PetscInt nrows,*rows; 108 PetscErrorCode ierr; 109 110 PetscFunctionBegin; 111 *zrows = NULL; 112 ierr = MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows);CHKERRQ(ierr); 113 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows);CHKERRQ(ierr); 114 PetscFunctionReturn(0); 115 } 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "MatFindNonzeroRows_SeqAIJ" 119 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 120 { 121 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 122 const MatScalar *aa; 123 PetscInt m=A->rmap->n,cnt = 0; 124 const PetscInt *ii; 125 PetscInt n,i,j,*rows; 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 *keptrows = 0; 130 ii = a->i; 131 for (i=0; i<m; i++) { 132 n = ii[i+1] - ii[i]; 133 if (!n) { 134 cnt++; 135 goto ok1; 136 } 137 aa = a->a + ii[i]; 138 for (j=0; j<n; j++) { 139 if (aa[j] != 0.0) goto ok1; 140 } 141 cnt++; 142 ok1:; 143 } 144 if (!cnt) PetscFunctionReturn(0); 145 ierr = PetscMalloc1(A->rmap->n-cnt,&rows);CHKERRQ(ierr); 146 cnt = 0; 147 for (i=0; i<m; i++) { 148 n = ii[i+1] - ii[i]; 149 if (!n) continue; 150 aa = a->a + ii[i]; 151 for (j=0; j<n; j++) { 152 if (aa[j] != 0.0) { 153 rows[cnt++] = i; 154 break; 155 } 156 } 157 } 158 ierr = ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows);CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 164 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 165 { 166 PetscErrorCode ierr; 167 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 168 PetscInt i,m = Y->rmap->n; 169 const PetscInt *diag; 170 MatScalar *aa = aij->a; 171 const PetscScalar *v; 172 PetscBool missing; 173 174 PetscFunctionBegin; 175 if (Y->assembled) { 176 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,NULL);CHKERRQ(ierr); 177 if (!missing) { 178 diag = aij->diag; 179 ierr = VecGetArrayRead(D,&v);CHKERRQ(ierr); 180 if (is == INSERT_VALUES) { 181 for (i=0; i<m; i++) { 182 aa[diag[i]] = v[i]; 183 } 184 } else { 185 for (i=0; i<m; i++) { 186 aa[diag[i]] += v[i]; 187 } 188 } 189 ierr = VecRestoreArrayRead(D,&v);CHKERRQ(ierr); 190 PetscFunctionReturn(0); 191 } 192 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 193 } 194 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 200 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 201 { 202 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 203 PetscErrorCode ierr; 204 PetscInt i,ishift; 205 206 PetscFunctionBegin; 207 *m = A->rmap->n; 208 if (!ia) PetscFunctionReturn(0); 209 ishift = 0; 210 if (symmetric && !A->structurally_symmetric) { 211 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 212 } else if (oshift == 1) { 213 PetscInt *tia; 214 PetscInt nz = a->i[A->rmap->n]; 215 /* malloc space and add 1 to i and j indices */ 216 ierr = PetscMalloc1(A->rmap->n+1,&tia);CHKERRQ(ierr); 217 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 218 *ia = tia; 219 if (ja) { 220 PetscInt *tja; 221 ierr = PetscMalloc1(nz+1,&tja);CHKERRQ(ierr); 222 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 223 *ja = tja; 224 } 225 } else { 226 *ia = a->i; 227 if (ja) *ja = a->j; 228 } 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 234 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 235 { 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 if (!ia) PetscFunctionReturn(0); 240 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 241 ierr = PetscFree(*ia);CHKERRQ(ierr); 242 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 243 } 244 PetscFunctionReturn(0); 245 } 246 247 #undef __FUNCT__ 248 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 249 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 250 { 251 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 252 PetscErrorCode ierr; 253 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 254 PetscInt nz = a->i[m],row,*jj,mr,col; 255 256 PetscFunctionBegin; 257 *nn = n; 258 if (!ia) PetscFunctionReturn(0); 259 if (symmetric) { 260 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,0,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 261 } else { 262 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 263 ierr = PetscMalloc1(n+1,&cia);CHKERRQ(ierr); 264 ierr = PetscMalloc1(nz+1,&cja);CHKERRQ(ierr); 265 jj = a->j; 266 for (i=0; i<nz; i++) { 267 collengths[jj[i]]++; 268 } 269 cia[0] = oshift; 270 for (i=0; i<n; i++) { 271 cia[i+1] = cia[i] + collengths[i]; 272 } 273 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 274 jj = a->j; 275 for (row=0; row<m; row++) { 276 mr = a->i[row+1] - a->i[row]; 277 for (i=0; i<mr; i++) { 278 col = *jj++; 279 280 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 281 } 282 } 283 ierr = PetscFree(collengths);CHKERRQ(ierr); 284 *ia = cia; *ja = cja; 285 } 286 PetscFunctionReturn(0); 287 } 288 289 #undef __FUNCT__ 290 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 291 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 292 { 293 PetscErrorCode ierr; 294 295 PetscFunctionBegin; 296 if (!ia) PetscFunctionReturn(0); 297 298 ierr = PetscFree(*ia);CHKERRQ(ierr); 299 ierr = PetscFree(*ja);CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 303 /* 304 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 305 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 306 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 307 */ 308 #undef __FUNCT__ 309 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ_Color" 310 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 311 { 312 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 313 PetscErrorCode ierr; 314 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 315 PetscInt nz = a->i[m],row,*jj,mr,col; 316 PetscInt *cspidx; 317 318 PetscFunctionBegin; 319 *nn = n; 320 if (!ia) PetscFunctionReturn(0); 321 322 ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr); 323 ierr = PetscMalloc1(n+1,&cia);CHKERRQ(ierr); 324 ierr = PetscMalloc1(nz+1,&cja);CHKERRQ(ierr); 325 ierr = PetscMalloc1(nz+1,&cspidx);CHKERRQ(ierr); 326 jj = a->j; 327 for (i=0; i<nz; i++) { 328 collengths[jj[i]]++; 329 } 330 cia[0] = oshift; 331 for (i=0; i<n; i++) { 332 cia[i+1] = cia[i] + collengths[i]; 333 } 334 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 335 jj = a->j; 336 for (row=0; row<m; row++) { 337 mr = a->i[row+1] - a->i[row]; 338 for (i=0; i<mr; i++) { 339 col = *jj++; 340 cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */ 341 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 342 } 343 } 344 ierr = PetscFree(collengths);CHKERRQ(ierr); 345 *ia = cia; *ja = cja; 346 *spidx = cspidx; 347 PetscFunctionReturn(0); 348 } 349 350 #undef __FUNCT__ 351 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ_Color" 352 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 353 { 354 PetscErrorCode ierr; 355 356 PetscFunctionBegin; 357 ierr = MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done);CHKERRQ(ierr); 358 ierr = PetscFree(*spidx);CHKERRQ(ierr); 359 PetscFunctionReturn(0); 360 } 361 362 #undef __FUNCT__ 363 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 364 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 365 { 366 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 367 PetscInt *ai = a->i; 368 PetscErrorCode ierr; 369 370 PetscFunctionBegin; 371 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 372 PetscFunctionReturn(0); 373 } 374 375 /* 376 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 377 378 - a single row of values is set with each call 379 - no row or column indices are negative or (in error) larger than the number of rows or columns 380 - the values are always added to the matrix, not set 381 - no new locations are introduced in the nonzero structure of the matrix 382 383 This does NOT assume the global column indices are sorted 384 385 */ 386 387 #include <petsc/private/isimpl.h> 388 #undef __FUNCT__ 389 #define __FUNCT__ "MatSeqAIJSetValuesLocalFast" 390 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 391 { 392 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 393 PetscInt low,high,t,row,nrow,i,col,l; 394 const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j; 395 PetscInt lastcol = -1; 396 MatScalar *ap,value,*aa = a->a; 397 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 398 399 row = ridx[im[0]]; 400 rp = aj + ai[row]; 401 ap = aa + ai[row]; 402 nrow = ailen[row]; 403 low = 0; 404 high = nrow; 405 for (l=0; l<n; l++) { /* loop over added columns */ 406 col = cidx[in[l]]; 407 value = v[l]; 408 409 if (col <= lastcol) low = 0; 410 else high = nrow; 411 lastcol = col; 412 while (high-low > 5) { 413 t = (low+high)/2; 414 if (rp[t] > col) high = t; 415 else low = t; 416 } 417 for (i=low; i<high; i++) { 418 if (rp[i] == col) { 419 ap[i] += value; 420 low = i + 1; 421 break; 422 } 423 } 424 } 425 return 0; 426 } 427 428 #undef __FUNCT__ 429 #define __FUNCT__ "MatSetValues_SeqAIJ" 430 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 431 { 432 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 433 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 434 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 435 PetscErrorCode ierr; 436 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 437 MatScalar *ap,value,*aa = a->a; 438 PetscBool ignorezeroentries = a->ignorezeroentries; 439 PetscBool roworiented = a->roworiented; 440 441 PetscFunctionBegin; 442 for (k=0; k<m; k++) { /* loop over added rows */ 443 row = im[k]; 444 if (row < 0) continue; 445 #if defined(PETSC_USE_DEBUG) 446 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); 447 #endif 448 rp = aj + ai[row]; ap = aa + ai[row]; 449 rmax = imax[row]; nrow = ailen[row]; 450 low = 0; 451 high = nrow; 452 for (l=0; l<n; l++) { /* loop over added columns */ 453 if (in[l] < 0) continue; 454 #if defined(PETSC_USE_DEBUG) 455 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); 456 #endif 457 col = in[l]; 458 if (roworiented) { 459 value = v[l + k*n]; 460 } else { 461 value = v[k + l*m]; 462 } 463 if ((value == 0.0 && ignorezeroentries) && (is == ADD_VALUES)) continue; 464 465 if (col <= lastcol) low = 0; 466 else high = nrow; 467 lastcol = col; 468 while (high-low > 5) { 469 t = (low+high)/2; 470 if (rp[t] > col) high = t; 471 else low = t; 472 } 473 for (i=low; i<high; i++) { 474 if (rp[i] > col) break; 475 if (rp[i] == col) { 476 if (is == ADD_VALUES) ap[i] += value; 477 else ap[i] = value; 478 low = i + 1; 479 goto noinsert; 480 } 481 } 482 if (value == 0.0 && ignorezeroentries) goto noinsert; 483 if (nonew == 1) goto noinsert; 484 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 485 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 486 N = nrow++ - 1; a->nz++; high++; 487 /* shift up all the later entries in this row */ 488 for (ii=N; ii>=i; ii--) { 489 rp[ii+1] = rp[ii]; 490 ap[ii+1] = ap[ii]; 491 } 492 rp[i] = col; 493 ap[i] = value; 494 low = i + 1; 495 A->nonzerostate++; 496 noinsert:; 497 } 498 ailen[row] = nrow; 499 } 500 PetscFunctionReturn(0); 501 } 502 503 504 #undef __FUNCT__ 505 #define __FUNCT__ "MatGetValues_SeqAIJ" 506 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 507 { 508 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 509 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 510 PetscInt *ai = a->i,*ailen = a->ilen; 511 MatScalar *ap,*aa = a->a; 512 513 PetscFunctionBegin; 514 for (k=0; k<m; k++) { /* loop over rows */ 515 row = im[k]; 516 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 517 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); 518 rp = aj + ai[row]; ap = aa + ai[row]; 519 nrow = ailen[row]; 520 for (l=0; l<n; l++) { /* loop over columns */ 521 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 522 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); 523 col = in[l]; 524 high = nrow; low = 0; /* assume unsorted */ 525 while (high-low > 5) { 526 t = (low+high)/2; 527 if (rp[t] > col) high = t; 528 else low = t; 529 } 530 for (i=low; i<high; i++) { 531 if (rp[i] > col) break; 532 if (rp[i] == col) { 533 *v++ = ap[i]; 534 goto finished; 535 } 536 } 537 *v++ = 0.0; 538 finished:; 539 } 540 } 541 PetscFunctionReturn(0); 542 } 543 544 545 #undef __FUNCT__ 546 #define __FUNCT__ "MatView_SeqAIJ_Binary" 547 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 548 { 549 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 550 PetscErrorCode ierr; 551 PetscInt i,*col_lens; 552 int fd; 553 FILE *file; 554 555 PetscFunctionBegin; 556 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 557 ierr = PetscMalloc1(4+A->rmap->n,&col_lens);CHKERRQ(ierr); 558 559 col_lens[0] = MAT_FILE_CLASSID; 560 col_lens[1] = A->rmap->n; 561 col_lens[2] = A->cmap->n; 562 col_lens[3] = a->nz; 563 564 /* store lengths of each row and write (including header) to file */ 565 for (i=0; i<A->rmap->n; i++) { 566 col_lens[4+i] = a->i[i+1] - a->i[i]; 567 } 568 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 569 ierr = PetscFree(col_lens);CHKERRQ(ierr); 570 571 /* store column indices (zero start index) */ 572 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 573 574 /* store nonzero values */ 575 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 576 577 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 578 if (file) { 579 fprintf(file,"-matload_block_size %d\n",(int)PetscAbs(A->rmap->bs)); 580 } 581 PetscFunctionReturn(0); 582 } 583 584 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 585 586 #undef __FUNCT__ 587 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 588 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 589 { 590 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 591 PetscErrorCode ierr; 592 PetscInt i,j,m = A->rmap->n; 593 const char *name; 594 PetscViewerFormat format; 595 596 PetscFunctionBegin; 597 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 598 if (format == PETSC_VIEWER_ASCII_MATLAB) { 599 PetscInt nofinalvalue = 0; 600 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) { 601 /* Need a dummy value to ensure the dimension of the matrix. */ 602 nofinalvalue = 1; 603 } 604 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 605 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 606 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 607 #if defined(PETSC_USE_COMPLEX) 608 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,4);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 609 #else 610 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 611 #endif 612 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 613 614 for (i=0; i<m; i++) { 615 for (j=a->i[i]; j<a->i[i+1]; j++) { 616 #if defined(PETSC_USE_COMPLEX) 617 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); 618 #else 619 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+1,(double)a->a[j]);CHKERRQ(ierr); 620 #endif 621 } 622 } 623 if (nofinalvalue) { 624 #if defined(PETSC_USE_COMPLEX) 625 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e %18.16e\n",m,A->cmap->n,0.,0.);CHKERRQ(ierr); 626 #else 627 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 628 #endif 629 } 630 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 631 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 632 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 633 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 634 PetscFunctionReturn(0); 635 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 636 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 637 for (i=0; i<m; i++) { 638 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 639 for (j=a->i[i]; j<a->i[i+1]; j++) { 640 #if defined(PETSC_USE_COMPLEX) 641 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 642 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 643 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 644 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 645 } else if (PetscRealPart(a->a[j]) != 0.0) { 646 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 647 } 648 #else 649 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr);} 650 #endif 651 } 652 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 653 } 654 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 655 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 656 PetscInt nzd=0,fshift=1,*sptr; 657 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 658 ierr = PetscMalloc1(m+1,&sptr);CHKERRQ(ierr); 659 for (i=0; i<m; i++) { 660 sptr[i] = nzd+1; 661 for (j=a->i[i]; j<a->i[i+1]; j++) { 662 if (a->j[j] >= i) { 663 #if defined(PETSC_USE_COMPLEX) 664 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 665 #else 666 if (a->a[j] != 0.0) nzd++; 667 #endif 668 } 669 } 670 } 671 sptr[m] = nzd+1; 672 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 673 for (i=0; i<m+1; i+=6) { 674 if (i+4<m) { 675 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); 676 } else if (i+3<m) { 677 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); 678 } else if (i+2<m) { 679 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr); 680 } else if (i+1<m) { 681 ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr); 682 } else if (i<m) { 683 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr); 684 } else { 685 ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr); 686 } 687 } 688 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 689 ierr = PetscFree(sptr);CHKERRQ(ierr); 690 for (i=0; i<m; i++) { 691 for (j=a->i[i]; j<a->i[i+1]; j++) { 692 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 693 } 694 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 695 } 696 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 697 for (i=0; i<m; i++) { 698 for (j=a->i[i]; j<a->i[i+1]; j++) { 699 if (a->j[j] >= i) { 700 #if defined(PETSC_USE_COMPLEX) 701 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 702 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 703 } 704 #else 705 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j]);CHKERRQ(ierr);} 706 #endif 707 } 708 } 709 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 710 } 711 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 712 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 713 PetscInt cnt = 0,jcnt; 714 PetscScalar value; 715 #if defined(PETSC_USE_COMPLEX) 716 PetscBool realonly = PETSC_TRUE; 717 718 for (i=0; i<a->i[m]; i++) { 719 if (PetscImaginaryPart(a->a[i]) != 0.0) { 720 realonly = PETSC_FALSE; 721 break; 722 } 723 } 724 #endif 725 726 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 727 for (i=0; i<m; i++) { 728 jcnt = 0; 729 for (j=0; j<A->cmap->n; j++) { 730 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 731 value = a->a[cnt++]; 732 jcnt++; 733 } else { 734 value = 0.0; 735 } 736 #if defined(PETSC_USE_COMPLEX) 737 if (realonly) { 738 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value));CHKERRQ(ierr); 739 } else { 740 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value));CHKERRQ(ierr); 741 } 742 #else 743 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value);CHKERRQ(ierr); 744 #endif 745 } 746 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 747 } 748 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 749 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 750 PetscInt fshift=1; 751 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 752 #if defined(PETSC_USE_COMPLEX) 753 ierr = PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n");CHKERRQ(ierr); 754 #else 755 ierr = PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n");CHKERRQ(ierr); 756 #endif 757 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 758 for (i=0; i<m; i++) { 759 for (j=a->i[i]; j<a->i[i+1]; j++) { 760 #if defined(PETSC_USE_COMPLEX) 761 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); 762 #else 763 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j]);CHKERRQ(ierr); 764 #endif 765 } 766 } 767 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 768 } else { 769 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 770 if (A->factortype) { 771 for (i=0; i<m; i++) { 772 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 773 /* L part */ 774 for (j=a->i[i]; j<a->i[i+1]; j++) { 775 #if defined(PETSC_USE_COMPLEX) 776 if (PetscImaginaryPart(a->a[j]) > 0.0) { 777 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 778 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 779 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));CHKERRQ(ierr); 780 } else { 781 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 782 } 783 #else 784 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 785 #endif 786 } 787 /* diagonal */ 788 j = a->diag[i]; 789 #if defined(PETSC_USE_COMPLEX) 790 if (PetscImaginaryPart(a->a[j]) > 0.0) { 791 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); 792 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 793 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); 794 } else { 795 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]));CHKERRQ(ierr); 796 } 797 #else 798 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)(1.0/a->a[j]));CHKERRQ(ierr); 799 #endif 800 801 /* U part */ 802 for (j=a->diag[i+1]+1; j<a->diag[i]; j++) { 803 #if defined(PETSC_USE_COMPLEX) 804 if (PetscImaginaryPart(a->a[j]) > 0.0) { 805 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 806 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 807 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])));CHKERRQ(ierr); 808 } else { 809 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 810 } 811 #else 812 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 813 #endif 814 } 815 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 816 } 817 } else { 818 for (i=0; i<m; i++) { 819 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 820 for (j=a->i[i]; j<a->i[i+1]; j++) { 821 #if defined(PETSC_USE_COMPLEX) 822 if (PetscImaginaryPart(a->a[j]) > 0.0) { 823 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 824 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 825 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 826 } else { 827 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)PetscRealPart(a->a[j]));CHKERRQ(ierr); 828 } 829 #else 830 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j],(double)a->a[j]);CHKERRQ(ierr); 831 #endif 832 } 833 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 834 } 835 } 836 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 837 } 838 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 839 PetscFunctionReturn(0); 840 } 841 842 #include <petscdraw.h> 843 #undef __FUNCT__ 844 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 845 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 846 { 847 Mat A = (Mat) Aa; 848 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 849 PetscErrorCode ierr; 850 PetscInt i,j,m = A->rmap->n; 851 int color; 852 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 853 PetscViewer viewer; 854 PetscViewerFormat format; 855 856 PetscFunctionBegin; 857 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 858 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 859 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 860 861 /* loop over matrix elements drawing boxes */ 862 863 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 864 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 865 /* Blue for negative, Cyan for zero and Red for positive */ 866 color = PETSC_DRAW_BLUE; 867 for (i=0; i<m; i++) { 868 y_l = m - i - 1.0; y_r = y_l + 1.0; 869 for (j=a->i[i]; j<a->i[i+1]; j++) { 870 x_l = a->j[j]; x_r = x_l + 1.0; 871 if (PetscRealPart(a->a[j]) >= 0.) continue; 872 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 873 } 874 } 875 color = PETSC_DRAW_CYAN; 876 for (i=0; i<m; i++) { 877 y_l = m - i - 1.0; y_r = y_l + 1.0; 878 for (j=a->i[i]; j<a->i[i+1]; j++) { 879 x_l = a->j[j]; x_r = x_l + 1.0; 880 if (a->a[j] != 0.) continue; 881 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 882 } 883 } 884 color = PETSC_DRAW_RED; 885 for (i=0; i<m; i++) { 886 y_l = m - i - 1.0; y_r = y_l + 1.0; 887 for (j=a->i[i]; j<a->i[i+1]; j++) { 888 x_l = a->j[j]; x_r = x_l + 1.0; 889 if (PetscRealPart(a->a[j]) <= 0.) continue; 890 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 891 } 892 } 893 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 894 } else { 895 /* use contour shading to indicate magnitude of values */ 896 /* first determine max of all nonzero values */ 897 PetscReal minv = 0.0, maxv = 0.0; 898 PetscInt nz = a->nz, count = 0; 899 PetscDraw popup; 900 901 for (i=0; i<nz; i++) { 902 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 903 } 904 if (minv >= maxv) maxv = minv + PETSC_SMALL; 905 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 906 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 907 908 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 909 for (i=0; i<m; i++) { 910 y_l = m - i - 1.0; 911 y_r = y_l + 1.0; 912 for (j=a->i[i]; j<a->i[i+1]; j++) { 913 x_l = a->j[j]; 914 x_r = x_l + 1.0; 915 color = PetscDrawRealToColor(PetscAbsScalar(a->a[count]),minv,maxv); 916 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 917 count++; 918 } 919 } 920 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 921 } 922 PetscFunctionReturn(0); 923 } 924 925 #include <petscdraw.h> 926 #undef __FUNCT__ 927 #define __FUNCT__ "MatView_SeqAIJ_Draw" 928 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 929 { 930 PetscErrorCode ierr; 931 PetscDraw draw; 932 PetscReal xr,yr,xl,yl,h,w; 933 PetscBool isnull; 934 935 PetscFunctionBegin; 936 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 937 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 938 if (isnull) PetscFunctionReturn(0); 939 940 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 941 xr += w; yr += h; xl = -w; yl = -h; 942 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 943 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 944 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 945 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 946 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 947 PetscFunctionReturn(0); 948 } 949 950 #undef __FUNCT__ 951 #define __FUNCT__ "MatView_SeqAIJ" 952 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 953 { 954 PetscErrorCode ierr; 955 PetscBool iascii,isbinary,isdraw; 956 957 PetscFunctionBegin; 958 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 959 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 960 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 961 if (iascii) { 962 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 963 } else if (isbinary) { 964 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 965 } else if (isdraw) { 966 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 967 } 968 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 969 PetscFunctionReturn(0); 970 } 971 972 #undef __FUNCT__ 973 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 974 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 975 { 976 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 977 PetscErrorCode ierr; 978 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 979 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 980 MatScalar *aa = a->a,*ap; 981 PetscReal ratio = 0.6; 982 983 PetscFunctionBegin; 984 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 985 986 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 987 for (i=1; i<m; i++) { 988 /* move each row back by the amount of empty slots (fshift) before it*/ 989 fshift += imax[i-1] - ailen[i-1]; 990 rmax = PetscMax(rmax,ailen[i]); 991 if (fshift) { 992 ip = aj + ai[i]; 993 ap = aa + ai[i]; 994 N = ailen[i]; 995 for (j=0; j<N; j++) { 996 ip[j-fshift] = ip[j]; 997 ap[j-fshift] = ap[j]; 998 } 999 } 1000 ai[i] = ai[i-1] + ailen[i-1]; 1001 } 1002 if (m) { 1003 fshift += imax[m-1] - ailen[m-1]; 1004 ai[m] = ai[m-1] + ailen[m-1]; 1005 } 1006 1007 /* reset ilen and imax for each row */ 1008 a->nonzerorowcnt = 0; 1009 for (i=0; i<m; i++) { 1010 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1011 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 1012 } 1013 a->nz = ai[m]; 1014 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); 1015 1016 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1017 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 1018 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 1019 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 1020 1021 A->info.mallocs += a->reallocs; 1022 a->reallocs = 0; 1023 A->info.nz_unneeded = (PetscReal)fshift; 1024 a->rmax = rmax; 1025 1026 ierr = MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 1027 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 1028 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1029 PetscFunctionReturn(0); 1030 } 1031 1032 #undef __FUNCT__ 1033 #define __FUNCT__ "MatRealPart_SeqAIJ" 1034 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1035 { 1036 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1037 PetscInt i,nz = a->nz; 1038 MatScalar *aa = a->a; 1039 PetscErrorCode ierr; 1040 1041 PetscFunctionBegin; 1042 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1043 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1044 PetscFunctionReturn(0); 1045 } 1046 1047 #undef __FUNCT__ 1048 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 1049 PetscErrorCode MatImaginaryPart_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] = PetscImaginaryPart(aa[i]); 1058 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1059 PetscFunctionReturn(0); 1060 } 1061 1062 #undef __FUNCT__ 1063 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 1064 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1065 { 1066 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1067 PetscErrorCode ierr; 1068 1069 PetscFunctionBegin; 1070 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 1071 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1072 PetscFunctionReturn(0); 1073 } 1074 1075 #undef __FUNCT__ 1076 #define __FUNCT__ "MatDestroy_SeqAIJ" 1077 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1078 { 1079 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1080 PetscErrorCode ierr; 1081 1082 PetscFunctionBegin; 1083 #if defined(PETSC_USE_LOG) 1084 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 1085 #endif 1086 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 1087 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1088 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1089 ierr = PetscFree(a->diag);CHKERRQ(ierr); 1090 ierr = PetscFree(a->ibdiag);CHKERRQ(ierr); 1091 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 1092 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 1093 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 1094 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 1095 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 1096 ierr = ISColoringDestroy(&a->coloring);CHKERRQ(ierr); 1097 ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);CHKERRQ(ierr); 1098 ierr = PetscFree(a->matmult_abdense);CHKERRQ(ierr); 1099 1100 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 1101 ierr = PetscFree(A->data);CHKERRQ(ierr); 1102 1103 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1104 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL);CHKERRQ(ierr); 1105 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);CHKERRQ(ierr); 1106 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);CHKERRQ(ierr); 1107 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL);CHKERRQ(ierr); 1108 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL);CHKERRQ(ierr); 1109 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL);CHKERRQ(ierr); 1110 #if defined(PETSC_HAVE_ELEMENTAL) 1111 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL);CHKERRQ(ierr); 1112 #endif 1113 #if defined(PETSC_HAVE_HYPRE) 1114 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_hypre_C",NULL);CHKERRQ(ierr); 1115 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMatMatMult_transpose_seqaij_seqaij_C",NULL);CHKERRQ(ierr); 1116 #endif 1117 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1118 ierr = PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);CHKERRQ(ierr); 1119 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);CHKERRQ(ierr); 1120 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr); 1121 ierr = PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);CHKERRQ(ierr); 1122 PetscFunctionReturn(0); 1123 } 1124 1125 #undef __FUNCT__ 1126 #define __FUNCT__ "MatSetOption_SeqAIJ" 1127 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1128 { 1129 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1130 PetscErrorCode ierr; 1131 1132 PetscFunctionBegin; 1133 switch (op) { 1134 case MAT_ROW_ORIENTED: 1135 a->roworiented = flg; 1136 break; 1137 case MAT_KEEP_NONZERO_PATTERN: 1138 a->keepnonzeropattern = flg; 1139 break; 1140 case MAT_NEW_NONZERO_LOCATIONS: 1141 a->nonew = (flg ? 0 : 1); 1142 break; 1143 case MAT_NEW_NONZERO_LOCATION_ERR: 1144 a->nonew = (flg ? -1 : 0); 1145 break; 1146 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1147 a->nonew = (flg ? -2 : 0); 1148 break; 1149 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1150 a->nounused = (flg ? -1 : 0); 1151 break; 1152 case MAT_IGNORE_ZERO_ENTRIES: 1153 a->ignorezeroentries = flg; 1154 break; 1155 case MAT_SPD: 1156 case MAT_SYMMETRIC: 1157 case MAT_STRUCTURALLY_SYMMETRIC: 1158 case MAT_HERMITIAN: 1159 case MAT_SYMMETRY_ETERNAL: 1160 /* These options are handled directly by MatSetOption() */ 1161 break; 1162 case MAT_NEW_DIAGONALS: 1163 case MAT_IGNORE_OFF_PROC_ENTRIES: 1164 case MAT_USE_HASH_TABLE: 1165 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1166 break; 1167 case MAT_USE_INODES: 1168 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 1169 break; 1170 case MAT_SUBMAT_SINGLEIS: 1171 A->submat_singleis = flg; 1172 break; 1173 default: 1174 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1175 } 1176 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 1177 PetscFunctionReturn(0); 1178 } 1179 1180 #undef __FUNCT__ 1181 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 1182 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1183 { 1184 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1185 PetscErrorCode ierr; 1186 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 1187 PetscScalar *aa=a->a,*x,zero=0.0; 1188 1189 PetscFunctionBegin; 1190 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1191 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1192 1193 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1194 PetscInt *diag=a->diag; 1195 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1196 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1197 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1198 PetscFunctionReturn(0); 1199 } 1200 1201 ierr = VecSet(v,zero);CHKERRQ(ierr); 1202 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1203 for (i=0; i<n; i++) { 1204 nz = ai[i+1] - ai[i]; 1205 if (!nz) x[i] = 0.0; 1206 for (j=ai[i]; j<ai[i+1]; j++) { 1207 if (aj[j] == i) { 1208 x[i] = aa[j]; 1209 break; 1210 } 1211 } 1212 } 1213 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1214 PetscFunctionReturn(0); 1215 } 1216 1217 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1218 #undef __FUNCT__ 1219 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 1220 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1221 { 1222 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1223 PetscScalar *y; 1224 const PetscScalar *x; 1225 PetscErrorCode ierr; 1226 PetscInt m = A->rmap->n; 1227 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1228 const MatScalar *v; 1229 PetscScalar alpha; 1230 PetscInt n,i,j; 1231 const PetscInt *idx,*ii,*ridx=NULL; 1232 Mat_CompressedRow cprow = a->compressedrow; 1233 PetscBool usecprow = cprow.use; 1234 #endif 1235 1236 PetscFunctionBegin; 1237 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 1238 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1239 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1240 1241 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1242 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 1243 #else 1244 if (usecprow) { 1245 m = cprow.nrows; 1246 ii = cprow.i; 1247 ridx = cprow.rindex; 1248 } else { 1249 ii = a->i; 1250 } 1251 for (i=0; i<m; i++) { 1252 idx = a->j + ii[i]; 1253 v = a->a + ii[i]; 1254 n = ii[i+1] - ii[i]; 1255 if (usecprow) { 1256 alpha = x[ridx[i]]; 1257 } else { 1258 alpha = x[i]; 1259 } 1260 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1261 } 1262 #endif 1263 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1264 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1265 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1266 PetscFunctionReturn(0); 1267 } 1268 1269 #undef __FUNCT__ 1270 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 1271 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1272 { 1273 PetscErrorCode ierr; 1274 1275 PetscFunctionBegin; 1276 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 1277 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 1278 PetscFunctionReturn(0); 1279 } 1280 1281 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1282 1283 #undef __FUNCT__ 1284 #define __FUNCT__ "MatMult_SeqAIJ" 1285 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1286 { 1287 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1288 PetscScalar *y; 1289 const PetscScalar *x; 1290 const MatScalar *aa; 1291 PetscErrorCode ierr; 1292 PetscInt m=A->rmap->n; 1293 const PetscInt *aj,*ii,*ridx=NULL; 1294 PetscInt n,i; 1295 PetscScalar sum; 1296 PetscBool usecprow=a->compressedrow.use; 1297 1298 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1299 #pragma disjoint(*x,*y,*aa) 1300 #endif 1301 1302 PetscFunctionBegin; 1303 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1304 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1305 ii = a->i; 1306 if (usecprow) { /* use compressed row format */ 1307 ierr = PetscMemzero(y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1308 m = a->compressedrow.nrows; 1309 ii = a->compressedrow.i; 1310 ridx = a->compressedrow.rindex; 1311 for (i=0; i<m; i++) { 1312 n = ii[i+1] - ii[i]; 1313 aj = a->j + ii[i]; 1314 aa = a->a + ii[i]; 1315 sum = 0.0; 1316 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1317 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1318 y[*ridx++] = sum; 1319 } 1320 } else { /* do not use compressed row format */ 1321 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1322 aj = a->j; 1323 aa = a->a; 1324 fortranmultaij_(&m,x,ii,aj,aa,y); 1325 #else 1326 for (i=0; i<m; i++) { 1327 n = ii[i+1] - ii[i]; 1328 aj = a->j + ii[i]; 1329 aa = a->a + ii[i]; 1330 sum = 0.0; 1331 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1332 y[i] = sum; 1333 } 1334 #endif 1335 } 1336 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1337 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1338 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1339 PetscFunctionReturn(0); 1340 } 1341 1342 #undef __FUNCT__ 1343 #define __FUNCT__ "MatMultMax_SeqAIJ" 1344 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1345 { 1346 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1347 PetscScalar *y; 1348 const PetscScalar *x; 1349 const MatScalar *aa; 1350 PetscErrorCode ierr; 1351 PetscInt m=A->rmap->n; 1352 const PetscInt *aj,*ii,*ridx=NULL; 1353 PetscInt n,i,nonzerorow=0; 1354 PetscScalar sum; 1355 PetscBool usecprow=a->compressedrow.use; 1356 1357 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1358 #pragma disjoint(*x,*y,*aa) 1359 #endif 1360 1361 PetscFunctionBegin; 1362 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1363 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1364 if (usecprow) { /* use compressed row format */ 1365 m = a->compressedrow.nrows; 1366 ii = a->compressedrow.i; 1367 ridx = a->compressedrow.rindex; 1368 for (i=0; i<m; i++) { 1369 n = ii[i+1] - ii[i]; 1370 aj = a->j + ii[i]; 1371 aa = a->a + ii[i]; 1372 sum = 0.0; 1373 nonzerorow += (n>0); 1374 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1375 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1376 y[*ridx++] = sum; 1377 } 1378 } else { /* do not use compressed row format */ 1379 ii = a->i; 1380 for (i=0; i<m; i++) { 1381 n = ii[i+1] - ii[i]; 1382 aj = a->j + ii[i]; 1383 aa = a->a + ii[i]; 1384 sum = 0.0; 1385 nonzerorow += (n>0); 1386 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1387 y[i] = sum; 1388 } 1389 } 1390 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1391 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1392 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1393 PetscFunctionReturn(0); 1394 } 1395 1396 #undef __FUNCT__ 1397 #define __FUNCT__ "MatMultAddMax_SeqAIJ" 1398 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1399 { 1400 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1401 PetscScalar *y,*z; 1402 const PetscScalar *x; 1403 const MatScalar *aa; 1404 PetscErrorCode ierr; 1405 PetscInt m = A->rmap->n,*aj,*ii; 1406 PetscInt n,i,*ridx=NULL; 1407 PetscScalar sum; 1408 PetscBool usecprow=a->compressedrow.use; 1409 1410 PetscFunctionBegin; 1411 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1412 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1413 if (usecprow) { /* use compressed row format */ 1414 if (zz != yy) { 1415 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1416 } 1417 m = a->compressedrow.nrows; 1418 ii = a->compressedrow.i; 1419 ridx = a->compressedrow.rindex; 1420 for (i=0; i<m; i++) { 1421 n = ii[i+1] - ii[i]; 1422 aj = a->j + ii[i]; 1423 aa = a->a + ii[i]; 1424 sum = y[*ridx]; 1425 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1426 z[*ridx++] = sum; 1427 } 1428 } else { /* do not use compressed row format */ 1429 ii = a->i; 1430 for (i=0; i<m; i++) { 1431 n = ii[i+1] - ii[i]; 1432 aj = a->j + ii[i]; 1433 aa = a->a + ii[i]; 1434 sum = y[i]; 1435 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1436 z[i] = sum; 1437 } 1438 } 1439 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1440 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1441 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1442 PetscFunctionReturn(0); 1443 } 1444 1445 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1446 #undef __FUNCT__ 1447 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1448 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1449 { 1450 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1451 PetscScalar *y,*z; 1452 const PetscScalar *x; 1453 const MatScalar *aa; 1454 PetscErrorCode ierr; 1455 const PetscInt *aj,*ii,*ridx=NULL; 1456 PetscInt m = A->rmap->n,n,i; 1457 PetscScalar sum; 1458 PetscBool usecprow=a->compressedrow.use; 1459 1460 PetscFunctionBegin; 1461 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1462 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1463 if (usecprow) { /* use compressed row format */ 1464 if (zz != yy) { 1465 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1466 } 1467 m = a->compressedrow.nrows; 1468 ii = a->compressedrow.i; 1469 ridx = a->compressedrow.rindex; 1470 for (i=0; i<m; i++) { 1471 n = ii[i+1] - ii[i]; 1472 aj = a->j + ii[i]; 1473 aa = a->a + ii[i]; 1474 sum = y[*ridx]; 1475 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1476 z[*ridx++] = sum; 1477 } 1478 } else { /* do not use compressed row format */ 1479 ii = a->i; 1480 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1481 aj = a->j; 1482 aa = a->a; 1483 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1484 #else 1485 for (i=0; i<m; i++) { 1486 n = ii[i+1] - ii[i]; 1487 aj = a->j + ii[i]; 1488 aa = a->a + ii[i]; 1489 sum = y[i]; 1490 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1491 z[i] = sum; 1492 } 1493 #endif 1494 } 1495 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1496 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1497 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1498 #if defined(PETSC_HAVE_CUSP) 1499 /* 1500 ierr = VecView(xx,0);CHKERRQ(ierr); 1501 ierr = VecView(zz,0);CHKERRQ(ierr); 1502 ierr = MatView(A,0);CHKERRQ(ierr); 1503 */ 1504 #endif 1505 PetscFunctionReturn(0); 1506 } 1507 1508 /* 1509 Adds diagonal pointers to sparse matrix structure. 1510 */ 1511 #undef __FUNCT__ 1512 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1513 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1514 { 1515 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1516 PetscErrorCode ierr; 1517 PetscInt i,j,m = A->rmap->n; 1518 1519 PetscFunctionBegin; 1520 if (!a->diag) { 1521 ierr = PetscMalloc1(m,&a->diag);CHKERRQ(ierr); 1522 ierr = PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));CHKERRQ(ierr); 1523 } 1524 for (i=0; i<A->rmap->n; i++) { 1525 a->diag[i] = a->i[i+1]; 1526 for (j=a->i[i]; j<a->i[i+1]; j++) { 1527 if (a->j[j] == i) { 1528 a->diag[i] = j; 1529 break; 1530 } 1531 } 1532 } 1533 PetscFunctionReturn(0); 1534 } 1535 1536 /* 1537 Checks for missing diagonals 1538 */ 1539 #undef __FUNCT__ 1540 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1541 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1542 { 1543 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1544 PetscInt *diag,*ii = a->i,i; 1545 1546 PetscFunctionBegin; 1547 *missing = PETSC_FALSE; 1548 if (A->rmap->n > 0 && !ii) { 1549 *missing = PETSC_TRUE; 1550 if (d) *d = 0; 1551 PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n"); 1552 } else { 1553 diag = a->diag; 1554 for (i=0; i<A->rmap->n; i++) { 1555 if (diag[i] >= ii[i+1]) { 1556 *missing = PETSC_TRUE; 1557 if (d) *d = i; 1558 PetscInfo1(A,"Matrix is missing diagonal number %D\n",i); 1559 break; 1560 } 1561 } 1562 } 1563 PetscFunctionReturn(0); 1564 } 1565 1566 #undef __FUNCT__ 1567 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1568 /* 1569 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1570 */ 1571 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1572 { 1573 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1574 PetscErrorCode ierr; 1575 PetscInt i,*diag,m = A->rmap->n; 1576 MatScalar *v = a->a; 1577 PetscScalar *idiag,*mdiag; 1578 1579 PetscFunctionBegin; 1580 if (a->idiagvalid) PetscFunctionReturn(0); 1581 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1582 diag = a->diag; 1583 if (!a->idiag) { 1584 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1585 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1586 v = a->a; 1587 } 1588 mdiag = a->mdiag; 1589 idiag = a->idiag; 1590 1591 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1592 for (i=0; i<m; i++) { 1593 mdiag[i] = v[diag[i]]; 1594 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1595 if (PetscRealPart(fshift)) { 1596 ierr = PetscInfo1(A,"Zero diagonal on row %D\n",i);CHKERRQ(ierr); 1597 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1598 A->factorerror_zeropivot_value = 0.0; 1599 A->factorerror_zeropivot_row = i; 1600 } SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1601 } 1602 idiag[i] = 1.0/v[diag[i]]; 1603 } 1604 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1605 } else { 1606 for (i=0; i<m; i++) { 1607 mdiag[i] = v[diag[i]]; 1608 idiag[i] = omega/(fshift + v[diag[i]]); 1609 } 1610 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1611 } 1612 a->idiagvalid = PETSC_TRUE; 1613 PetscFunctionReturn(0); 1614 } 1615 1616 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "MatSOR_SeqAIJ" 1619 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1620 { 1621 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1622 PetscScalar *x,d,sum,*t,scale; 1623 const MatScalar *v,*idiag=0,*mdiag; 1624 const PetscScalar *b, *bs,*xb, *ts; 1625 PetscErrorCode ierr; 1626 PetscInt n,m = A->rmap->n,i; 1627 const PetscInt *idx,*diag; 1628 1629 PetscFunctionBegin; 1630 its = its*lits; 1631 1632 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1633 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1634 a->fshift = fshift; 1635 a->omega = omega; 1636 1637 diag = a->diag; 1638 t = a->ssor_work; 1639 idiag = a->idiag; 1640 mdiag = a->mdiag; 1641 1642 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1643 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1644 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1645 if (flag == SOR_APPLY_UPPER) { 1646 /* apply (U + D/omega) to the vector */ 1647 bs = b; 1648 for (i=0; i<m; i++) { 1649 d = fshift + mdiag[i]; 1650 n = a->i[i+1] - diag[i] - 1; 1651 idx = a->j + diag[i] + 1; 1652 v = a->a + diag[i] + 1; 1653 sum = b[i]*d/omega; 1654 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1655 x[i] = sum; 1656 } 1657 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1658 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1659 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1660 PetscFunctionReturn(0); 1661 } 1662 1663 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1664 else if (flag & SOR_EISENSTAT) { 1665 /* Let A = L + U + D; where L is lower trianglar, 1666 U is upper triangular, E = D/omega; This routine applies 1667 1668 (L + E)^{-1} A (U + E)^{-1} 1669 1670 to a vector efficiently using Eisenstat's trick. 1671 */ 1672 scale = (2.0/omega) - 1.0; 1673 1674 /* x = (E + U)^{-1} b */ 1675 for (i=m-1; i>=0; i--) { 1676 n = a->i[i+1] - diag[i] - 1; 1677 idx = a->j + diag[i] + 1; 1678 v = a->a + diag[i] + 1; 1679 sum = b[i]; 1680 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1681 x[i] = sum*idiag[i]; 1682 } 1683 1684 /* t = b - (2*E - D)x */ 1685 v = a->a; 1686 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1687 1688 /* t = (E + L)^{-1}t */ 1689 ts = t; 1690 diag = a->diag; 1691 for (i=0; i<m; i++) { 1692 n = diag[i] - a->i[i]; 1693 idx = a->j + a->i[i]; 1694 v = a->a + a->i[i]; 1695 sum = t[i]; 1696 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1697 t[i] = sum*idiag[i]; 1698 /* x = x + t */ 1699 x[i] += t[i]; 1700 } 1701 1702 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1703 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1704 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1705 PetscFunctionReturn(0); 1706 } 1707 if (flag & SOR_ZERO_INITIAL_GUESS) { 1708 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1709 for (i=0; i<m; i++) { 1710 n = diag[i] - a->i[i]; 1711 idx = a->j + a->i[i]; 1712 v = a->a + a->i[i]; 1713 sum = b[i]; 1714 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1715 t[i] = sum; 1716 x[i] = sum*idiag[i]; 1717 } 1718 xb = t; 1719 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1720 } else xb = b; 1721 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1722 for (i=m-1; i>=0; i--) { 1723 n = a->i[i+1] - diag[i] - 1; 1724 idx = a->j + diag[i] + 1; 1725 v = a->a + diag[i] + 1; 1726 sum = xb[i]; 1727 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1728 if (xb == b) { 1729 x[i] = sum*idiag[i]; 1730 } else { 1731 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1732 } 1733 } 1734 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1735 } 1736 its--; 1737 } 1738 while (its--) { 1739 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1740 for (i=0; i<m; i++) { 1741 /* lower */ 1742 n = diag[i] - a->i[i]; 1743 idx = a->j + a->i[i]; 1744 v = a->a + a->i[i]; 1745 sum = b[i]; 1746 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1747 t[i] = sum; /* save application of the lower-triangular part */ 1748 /* upper */ 1749 n = a->i[i+1] - diag[i] - 1; 1750 idx = a->j + diag[i] + 1; 1751 v = a->a + diag[i] + 1; 1752 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1753 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1754 } 1755 xb = t; 1756 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1757 } else xb = b; 1758 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1759 for (i=m-1; i>=0; i--) { 1760 sum = xb[i]; 1761 if (xb == b) { 1762 /* whole matrix (no checkpointing available) */ 1763 n = a->i[i+1] - a->i[i]; 1764 idx = a->j + a->i[i]; 1765 v = a->a + a->i[i]; 1766 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1767 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1768 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1769 n = a->i[i+1] - diag[i] - 1; 1770 idx = a->j + diag[i] + 1; 1771 v = a->a + diag[i] + 1; 1772 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1773 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1774 } 1775 } 1776 if (xb == b) { 1777 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1778 } else { 1779 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1780 } 1781 } 1782 } 1783 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1784 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1785 PetscFunctionReturn(0); 1786 } 1787 1788 1789 #undef __FUNCT__ 1790 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1791 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1792 { 1793 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1794 1795 PetscFunctionBegin; 1796 info->block_size = 1.0; 1797 info->nz_allocated = (double)a->maxnz; 1798 info->nz_used = (double)a->nz; 1799 info->nz_unneeded = (double)(a->maxnz - a->nz); 1800 info->assemblies = (double)A->num_ass; 1801 info->mallocs = (double)A->info.mallocs; 1802 info->memory = ((PetscObject)A)->mem; 1803 if (A->factortype) { 1804 info->fill_ratio_given = A->info.fill_ratio_given; 1805 info->fill_ratio_needed = A->info.fill_ratio_needed; 1806 info->factor_mallocs = A->info.factor_mallocs; 1807 } else { 1808 info->fill_ratio_given = 0; 1809 info->fill_ratio_needed = 0; 1810 info->factor_mallocs = 0; 1811 } 1812 PetscFunctionReturn(0); 1813 } 1814 1815 #undef __FUNCT__ 1816 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1817 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1818 { 1819 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1820 PetscInt i,m = A->rmap->n - 1,d = 0; 1821 PetscErrorCode ierr; 1822 const PetscScalar *xx; 1823 PetscScalar *bb; 1824 PetscBool missing; 1825 1826 PetscFunctionBegin; 1827 if (x && b) { 1828 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1829 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1830 for (i=0; i<N; i++) { 1831 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1832 bb[rows[i]] = diag*xx[rows[i]]; 1833 } 1834 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1835 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1836 } 1837 1838 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1839 if (a->keepnonzeropattern) { 1840 for (i=0; i<N; i++) { 1841 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1842 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1843 } 1844 if (diag != 0.0) { 1845 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1846 for (i=0; i<N; i++) { 1847 a->a[a->diag[rows[i]]] = diag; 1848 } 1849 } 1850 } else { 1851 if (diag != 0.0) { 1852 for (i=0; i<N; i++) { 1853 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1854 if (a->ilen[rows[i]] > 0) { 1855 a->ilen[rows[i]] = 1; 1856 a->a[a->i[rows[i]]] = diag; 1857 a->j[a->i[rows[i]]] = rows[i]; 1858 } else { /* in case row was completely empty */ 1859 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1860 } 1861 } 1862 } else { 1863 for (i=0; i<N; i++) { 1864 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1865 a->ilen[rows[i]] = 0; 1866 } 1867 } 1868 A->nonzerostate++; 1869 } 1870 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1871 PetscFunctionReturn(0); 1872 } 1873 1874 #undef __FUNCT__ 1875 #define __FUNCT__ "MatZeroRowsColumns_SeqAIJ" 1876 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1877 { 1878 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1879 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1880 PetscErrorCode ierr; 1881 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 1882 const PetscScalar *xx; 1883 PetscScalar *bb; 1884 1885 PetscFunctionBegin; 1886 if (x && b) { 1887 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1888 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1889 vecs = PETSC_TRUE; 1890 } 1891 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1892 for (i=0; i<N; i++) { 1893 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1894 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1895 1896 zeroed[rows[i]] = PETSC_TRUE; 1897 } 1898 for (i=0; i<A->rmap->n; i++) { 1899 if (!zeroed[i]) { 1900 for (j=a->i[i]; j<a->i[i+1]; j++) { 1901 if (zeroed[a->j[j]]) { 1902 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 1903 a->a[j] = 0.0; 1904 } 1905 } 1906 } else if (vecs) bb[i] = diag*xx[i]; 1907 } 1908 if (x && b) { 1909 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1910 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1911 } 1912 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1913 if (diag != 0.0) { 1914 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1915 if (missing) { 1916 if (a->nonew) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1917 else { 1918 for (i=0; i<N; i++) { 1919 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1920 } 1921 } 1922 } else { 1923 for (i=0; i<N; i++) { 1924 a->a[a->diag[rows[i]]] = diag; 1925 } 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 0, 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 only 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 B->was_assembled = PETSC_FALSE; 3656 B->assembled = PETSC_FALSE; 3657 PetscFunctionReturn(0); 3658 } 3659 3660 #undef __FUNCT__ 3661 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3662 /*@ 3663 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3664 3665 Input Parameters: 3666 + B - the matrix 3667 . i - the indices into j for the start of each row (starts with zero) 3668 . j - the column indices for each row (starts with zero) these must be sorted for each row 3669 - v - optional values in the matrix 3670 3671 Level: developer 3672 3673 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3674 3675 .keywords: matrix, aij, compressed row, sparse, sequential 3676 3677 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3678 @*/ 3679 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3680 { 3681 PetscErrorCode ierr; 3682 3683 PetscFunctionBegin; 3684 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3685 PetscValidType(B,1); 3686 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3687 PetscFunctionReturn(0); 3688 } 3689 3690 #undef __FUNCT__ 3691 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3692 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3693 { 3694 PetscInt i; 3695 PetscInt m,n; 3696 PetscInt nz; 3697 PetscInt *nnz, nz_max = 0; 3698 PetscScalar *values; 3699 PetscErrorCode ierr; 3700 3701 PetscFunctionBegin; 3702 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3703 3704 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3705 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3706 3707 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3708 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3709 for (i = 0; i < m; i++) { 3710 nz = Ii[i+1]- Ii[i]; 3711 nz_max = PetscMax(nz_max, nz); 3712 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3713 nnz[i] = nz; 3714 } 3715 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3716 ierr = PetscFree(nnz);CHKERRQ(ierr); 3717 3718 if (v) { 3719 values = (PetscScalar*) v; 3720 } else { 3721 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3722 } 3723 3724 for (i = 0; i < m; i++) { 3725 nz = Ii[i+1] - Ii[i]; 3726 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3727 } 3728 3729 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3730 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3731 3732 if (!v) { 3733 ierr = PetscFree(values);CHKERRQ(ierr); 3734 } 3735 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3736 PetscFunctionReturn(0); 3737 } 3738 3739 #include <../src/mat/impls/dense/seq/dense.h> 3740 #include <petsc/private/kernels/petscaxpy.h> 3741 3742 #undef __FUNCT__ 3743 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3744 /* 3745 Computes (B'*A')' since computing B*A directly is untenable 3746 3747 n p p 3748 ( ) ( ) ( ) 3749 m ( A ) * n ( B ) = m ( C ) 3750 ( ) ( ) ( ) 3751 3752 */ 3753 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3754 { 3755 PetscErrorCode ierr; 3756 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3757 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3758 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3759 PetscInt i,n,m,q,p; 3760 const PetscInt *ii,*idx; 3761 const PetscScalar *b,*a,*a_q; 3762 PetscScalar *c,*c_q; 3763 3764 PetscFunctionBegin; 3765 m = A->rmap->n; 3766 n = A->cmap->n; 3767 p = B->cmap->n; 3768 a = sub_a->v; 3769 b = sub_b->a; 3770 c = sub_c->v; 3771 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3772 3773 ii = sub_b->i; 3774 idx = sub_b->j; 3775 for (i=0; i<n; i++) { 3776 q = ii[i+1] - ii[i]; 3777 while (q-->0) { 3778 c_q = c + m*(*idx); 3779 a_q = a + m*i; 3780 PetscKernelAXPY(c_q,*b,a_q,m); 3781 idx++; 3782 b++; 3783 } 3784 } 3785 PetscFunctionReturn(0); 3786 } 3787 3788 #undef __FUNCT__ 3789 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3790 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3791 { 3792 PetscErrorCode ierr; 3793 PetscInt m=A->rmap->n,n=B->cmap->n; 3794 Mat Cmat; 3795 3796 PetscFunctionBegin; 3797 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); 3798 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3799 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3800 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3801 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3802 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3803 3804 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3805 3806 *C = Cmat; 3807 PetscFunctionReturn(0); 3808 } 3809 3810 /* ----------------------------------------------------------------*/ 3811 #undef __FUNCT__ 3812 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3813 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3814 { 3815 PetscErrorCode ierr; 3816 3817 PetscFunctionBegin; 3818 if (scall == MAT_INITIAL_MATRIX) { 3819 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3820 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3821 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3822 } 3823 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3824 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3825 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3826 PetscFunctionReturn(0); 3827 } 3828 3829 3830 /*MC 3831 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3832 based on compressed sparse row format. 3833 3834 Options Database Keys: 3835 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3836 3837 Level: beginner 3838 3839 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3840 M*/ 3841 3842 /*MC 3843 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3844 3845 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3846 and MATMPIAIJ otherwise. As a result, for single process communicators, 3847 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3848 for communicators controlling multiple processes. It is recommended that you call both of 3849 the above preallocation routines for simplicity. 3850 3851 Options Database Keys: 3852 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3853 3854 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3855 enough exist. 3856 3857 Level: beginner 3858 3859 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3860 M*/ 3861 3862 /*MC 3863 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3864 3865 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3866 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3867 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3868 for communicators controlling multiple processes. It is recommended that you call both of 3869 the above preallocation routines for simplicity. 3870 3871 Options Database Keys: 3872 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3873 3874 Level: beginner 3875 3876 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3877 M*/ 3878 3879 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3880 #if defined(PETSC_HAVE_ELEMENTAL) 3881 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3882 #endif 3883 #if defined(PETSC_HAVE_HYPRE) 3884 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 3885 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 3886 #endif 3887 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 3888 3889 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3890 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3891 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3892 #endif 3893 3894 3895 #undef __FUNCT__ 3896 #define __FUNCT__ "MatSeqAIJGetArray" 3897 /*@C 3898 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 3899 3900 Not Collective 3901 3902 Input Parameter: 3903 . mat - a MATSEQAIJ matrix 3904 3905 Output Parameter: 3906 . array - pointer to the data 3907 3908 Level: intermediate 3909 3910 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3911 @*/ 3912 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3913 { 3914 PetscErrorCode ierr; 3915 3916 PetscFunctionBegin; 3917 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3918 PetscFunctionReturn(0); 3919 } 3920 3921 #undef __FUNCT__ 3922 #define __FUNCT__ "MatSeqAIJGetMaxRowNonzeros" 3923 /*@C 3924 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 3925 3926 Not Collective 3927 3928 Input Parameter: 3929 . mat - a MATSEQAIJ matrix 3930 3931 Output Parameter: 3932 . nz - the maximum number of nonzeros in any row 3933 3934 Level: intermediate 3935 3936 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3937 @*/ 3938 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 3939 { 3940 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 3941 3942 PetscFunctionBegin; 3943 *nz = aij->rmax; 3944 PetscFunctionReturn(0); 3945 } 3946 3947 #undef __FUNCT__ 3948 #define __FUNCT__ "MatSeqAIJRestoreArray" 3949 /*@C 3950 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 3951 3952 Not Collective 3953 3954 Input Parameters: 3955 . mat - a MATSEQAIJ matrix 3956 . array - pointer to the data 3957 3958 Level: intermediate 3959 3960 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 3961 @*/ 3962 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3963 { 3964 PetscErrorCode ierr; 3965 3966 PetscFunctionBegin; 3967 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3968 PetscFunctionReturn(0); 3969 } 3970 3971 #undef __FUNCT__ 3972 #define __FUNCT__ "MatCreate_SeqAIJ" 3973 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 3974 { 3975 Mat_SeqAIJ *b; 3976 PetscErrorCode ierr; 3977 PetscMPIInt size; 3978 3979 PetscFunctionBegin; 3980 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3981 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3982 3983 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3984 3985 B->data = (void*)b; 3986 3987 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3988 3989 b->row = 0; 3990 b->col = 0; 3991 b->icol = 0; 3992 b->reallocs = 0; 3993 b->ignorezeroentries = PETSC_FALSE; 3994 b->roworiented = PETSC_TRUE; 3995 b->nonew = 0; 3996 b->diag = 0; 3997 b->solve_work = 0; 3998 B->spptr = 0; 3999 b->saved_values = 0; 4000 b->idiag = 0; 4001 b->mdiag = 0; 4002 b->ssor_work = 0; 4003 b->omega = 1.0; 4004 b->fshift = 0.0; 4005 b->idiagvalid = PETSC_FALSE; 4006 b->ibdiagvalid = PETSC_FALSE; 4007 b->keepnonzeropattern = PETSC_FALSE; 4008 4009 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4010 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4011 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4012 4013 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4014 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4015 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4016 #endif 4017 4018 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4019 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4020 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4021 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4022 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4023 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4024 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4025 #if defined(PETSC_HAVE_ELEMENTAL) 4026 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4027 #endif 4028 #if defined(PETSC_HAVE_HYPRE) 4029 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 4030 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 4031 #endif 4032 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 4033 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4034 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4035 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4036 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4037 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4038 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4039 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4040 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4041 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4042 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4043 PetscFunctionReturn(0); 4044 } 4045 4046 #undef __FUNCT__ 4047 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 4048 /* 4049 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4050 */ 4051 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4052 { 4053 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4054 PetscErrorCode ierr; 4055 PetscInt i,m = A->rmap->n; 4056 4057 PetscFunctionBegin; 4058 c = (Mat_SeqAIJ*)C->data; 4059 4060 C->factortype = A->factortype; 4061 c->row = 0; 4062 c->col = 0; 4063 c->icol = 0; 4064 c->reallocs = 0; 4065 4066 C->assembled = PETSC_TRUE; 4067 4068 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4069 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4070 4071 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4072 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4073 for (i=0; i<m; i++) { 4074 c->imax[i] = a->imax[i]; 4075 c->ilen[i] = a->ilen[i]; 4076 } 4077 4078 /* allocate the matrix space */ 4079 if (mallocmatspace) { 4080 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4081 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4082 4083 c->singlemalloc = PETSC_TRUE; 4084 4085 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4086 if (m > 0) { 4087 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4088 if (cpvalues == MAT_COPY_VALUES) { 4089 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4090 } else { 4091 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4092 } 4093 } 4094 } 4095 4096 c->ignorezeroentries = a->ignorezeroentries; 4097 c->roworiented = a->roworiented; 4098 c->nonew = a->nonew; 4099 if (a->diag) { 4100 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4101 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4102 for (i=0; i<m; i++) { 4103 c->diag[i] = a->diag[i]; 4104 } 4105 } else c->diag = 0; 4106 4107 c->solve_work = 0; 4108 c->saved_values = 0; 4109 c->idiag = 0; 4110 c->ssor_work = 0; 4111 c->keepnonzeropattern = a->keepnonzeropattern; 4112 c->free_a = PETSC_TRUE; 4113 c->free_ij = PETSC_TRUE; 4114 4115 c->rmax = a->rmax; 4116 c->nz = a->nz; 4117 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4118 C->preallocated = PETSC_TRUE; 4119 4120 c->compressedrow.use = a->compressedrow.use; 4121 c->compressedrow.nrows = a->compressedrow.nrows; 4122 if (a->compressedrow.use) { 4123 i = a->compressedrow.nrows; 4124 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4125 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4126 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4127 } else { 4128 c->compressedrow.use = PETSC_FALSE; 4129 c->compressedrow.i = NULL; 4130 c->compressedrow.rindex = NULL; 4131 } 4132 c->nonzerorowcnt = a->nonzerorowcnt; 4133 C->nonzerostate = A->nonzerostate; 4134 4135 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4136 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4137 PetscFunctionReturn(0); 4138 } 4139 4140 #undef __FUNCT__ 4141 #define __FUNCT__ "MatDuplicate_SeqAIJ" 4142 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4143 { 4144 PetscErrorCode ierr; 4145 4146 PetscFunctionBegin; 4147 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4148 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4149 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4150 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4151 } 4152 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4153 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4154 PetscFunctionReturn(0); 4155 } 4156 4157 #undef __FUNCT__ 4158 #define __FUNCT__ "MatLoad_SeqAIJ" 4159 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4160 { 4161 Mat_SeqAIJ *a; 4162 PetscErrorCode ierr; 4163 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4164 int fd; 4165 PetscMPIInt size; 4166 MPI_Comm comm; 4167 PetscInt bs = newMat->rmap->bs; 4168 4169 PetscFunctionBegin; 4170 /* force binary viewer to load .info file if it has not yet done so */ 4171 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4172 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4173 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4174 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4175 4176 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4177 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4178 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4179 if (bs < 0) bs = 1; 4180 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4181 4182 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4183 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4184 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4185 M = header[1]; N = header[2]; nz = header[3]; 4186 4187 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4188 4189 /* read in row lengths */ 4190 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4191 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4192 4193 /* check if sum of rowlengths is same as nz */ 4194 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4195 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); 4196 4197 /* set global size if not set already*/ 4198 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4199 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4200 } else { 4201 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4202 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4203 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4204 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4205 } 4206 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); 4207 } 4208 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4209 a = (Mat_SeqAIJ*)newMat->data; 4210 4211 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4212 4213 /* read in nonzero values */ 4214 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4215 4216 /* set matrix "i" values */ 4217 a->i[0] = 0; 4218 for (i=1; i<= M; i++) { 4219 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4220 a->ilen[i-1] = rowlengths[i-1]; 4221 } 4222 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4223 4224 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4225 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4226 PetscFunctionReturn(0); 4227 } 4228 4229 #undef __FUNCT__ 4230 #define __FUNCT__ "MatEqual_SeqAIJ" 4231 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4232 { 4233 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4234 PetscErrorCode ierr; 4235 #if defined(PETSC_USE_COMPLEX) 4236 PetscInt k; 4237 #endif 4238 4239 PetscFunctionBegin; 4240 /* If the matrix dimensions are not equal,or no of nonzeros */ 4241 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4242 *flg = PETSC_FALSE; 4243 PetscFunctionReturn(0); 4244 } 4245 4246 /* if the a->i are the same */ 4247 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4248 if (!*flg) PetscFunctionReturn(0); 4249 4250 /* if a->j are the same */ 4251 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4252 if (!*flg) PetscFunctionReturn(0); 4253 4254 /* if a->a are the same */ 4255 #if defined(PETSC_USE_COMPLEX) 4256 for (k=0; k<a->nz; k++) { 4257 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4258 *flg = PETSC_FALSE; 4259 PetscFunctionReturn(0); 4260 } 4261 } 4262 #else 4263 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4264 #endif 4265 PetscFunctionReturn(0); 4266 } 4267 4268 #undef __FUNCT__ 4269 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4270 /*@ 4271 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4272 provided by the user. 4273 4274 Collective on MPI_Comm 4275 4276 Input Parameters: 4277 + comm - must be an MPI communicator of size 1 4278 . m - number of rows 4279 . n - number of columns 4280 . i - row indices 4281 . j - column indices 4282 - a - matrix values 4283 4284 Output Parameter: 4285 . mat - the matrix 4286 4287 Level: intermediate 4288 4289 Notes: 4290 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4291 once the matrix is destroyed and not before 4292 4293 You cannot set new nonzero locations into this matrix, that will generate an error. 4294 4295 The i and j indices are 0 based 4296 4297 The format which is used for the sparse matrix input, is equivalent to a 4298 row-major ordering.. i.e for the following matrix, the input data expected is 4299 as shown 4300 4301 $ 1 0 0 4302 $ 2 0 3 4303 $ 4 5 6 4304 $ 4305 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4306 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4307 $ v = {1,2,3,4,5,6} [size = 6] 4308 4309 4310 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4311 4312 @*/ 4313 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 4314 { 4315 PetscErrorCode ierr; 4316 PetscInt ii; 4317 Mat_SeqAIJ *aij; 4318 #if defined(PETSC_USE_DEBUG) 4319 PetscInt jj; 4320 #endif 4321 4322 PetscFunctionBegin; 4323 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4324 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4325 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4326 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4327 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4328 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4329 aij = (Mat_SeqAIJ*)(*mat)->data; 4330 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4331 4332 aij->i = i; 4333 aij->j = j; 4334 aij->a = a; 4335 aij->singlemalloc = PETSC_FALSE; 4336 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4337 aij->free_a = PETSC_FALSE; 4338 aij->free_ij = PETSC_FALSE; 4339 4340 for (ii=0; ii<m; ii++) { 4341 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4342 #if defined(PETSC_USE_DEBUG) 4343 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]); 4344 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4345 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); 4346 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); 4347 } 4348 #endif 4349 } 4350 #if defined(PETSC_USE_DEBUG) 4351 for (ii=0; ii<aij->i[m]; ii++) { 4352 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4353 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]); 4354 } 4355 #endif 4356 4357 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4358 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4359 PetscFunctionReturn(0); 4360 } 4361 #undef __FUNCT__ 4362 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4363 /*@C 4364 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4365 provided by the user. 4366 4367 Collective on MPI_Comm 4368 4369 Input Parameters: 4370 + comm - must be an MPI communicator of size 1 4371 . m - number of rows 4372 . n - number of columns 4373 . i - row indices 4374 . j - column indices 4375 . a - matrix values 4376 . nz - number of nonzeros 4377 - idx - 0 or 1 based 4378 4379 Output Parameter: 4380 . mat - the matrix 4381 4382 Level: intermediate 4383 4384 Notes: 4385 The i and j indices are 0 based 4386 4387 The format which is used for the sparse matrix input, is equivalent to a 4388 row-major ordering.. i.e for the following matrix, the input data expected is 4389 as shown: 4390 4391 1 0 0 4392 2 0 3 4393 4 5 6 4394 4395 i = {0,1,1,2,2,2} 4396 j = {0,0,2,0,1,2} 4397 v = {1,2,3,4,5,6} 4398 4399 4400 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4401 4402 @*/ 4403 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4404 { 4405 PetscErrorCode ierr; 4406 PetscInt ii, *nnz, one = 1,row,col; 4407 4408 4409 PetscFunctionBegin; 4410 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4411 for (ii = 0; ii < nz; ii++) { 4412 nnz[i[ii] - !!idx] += 1; 4413 } 4414 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4415 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4416 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4417 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4418 for (ii = 0; ii < nz; ii++) { 4419 if (idx) { 4420 row = i[ii] - 1; 4421 col = j[ii] - 1; 4422 } else { 4423 row = i[ii]; 4424 col = j[ii]; 4425 } 4426 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4427 } 4428 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4429 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4430 ierr = PetscFree(nnz);CHKERRQ(ierr); 4431 PetscFunctionReturn(0); 4432 } 4433 4434 #undef __FUNCT__ 4435 #define __FUNCT__ "MatSeqAIJInvalidateDiagonal" 4436 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4437 { 4438 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4439 PetscErrorCode ierr; 4440 4441 PetscFunctionBegin; 4442 a->idiagvalid = PETSC_FALSE; 4443 a->ibdiagvalid = PETSC_FALSE; 4444 4445 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4446 PetscFunctionReturn(0); 4447 } 4448 4449 #undef __FUNCT__ 4450 #define __FUNCT__ "MatCreateMPIMatConcatenateSeqMat_SeqAIJ" 4451 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4452 { 4453 PetscErrorCode ierr; 4454 4455 PetscFunctionBegin; 4456 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4457 PetscFunctionReturn(0); 4458 } 4459 4460 /* 4461 Permute A into C's *local* index space using rowemb,colemb. 4462 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4463 of [0,m), colemb is in [0,n). 4464 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4465 */ 4466 #undef __FUNCT__ 4467 #define __FUNCT__ "MatSetSeqMat_SeqAIJ" 4468 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4469 { 4470 /* If making this function public, change the error returned in this function away from _PLIB. */ 4471 PetscErrorCode ierr; 4472 Mat_SeqAIJ *Baij; 4473 PetscBool seqaij; 4474 PetscInt m,n,*nz,i,j,count; 4475 PetscScalar v; 4476 const PetscInt *rowindices,*colindices; 4477 4478 PetscFunctionBegin; 4479 if (!B) PetscFunctionReturn(0); 4480 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4481 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4482 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4483 if (rowemb) { 4484 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4485 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); 4486 } else { 4487 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4488 } 4489 if (colemb) { 4490 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4491 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); 4492 } else { 4493 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4494 } 4495 4496 Baij = (Mat_SeqAIJ*)(B->data); 4497 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4498 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4499 for (i=0; i<B->rmap->n; i++) { 4500 nz[i] = Baij->i[i+1] - Baij->i[i]; 4501 } 4502 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4503 ierr = PetscFree(nz);CHKERRQ(ierr); 4504 } 4505 if (pattern == SUBSET_NONZERO_PATTERN) { 4506 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4507 } 4508 count = 0; 4509 rowindices = NULL; 4510 colindices = NULL; 4511 if (rowemb) { 4512 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4513 } 4514 if (colemb) { 4515 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4516 } 4517 for (i=0; i<B->rmap->n; i++) { 4518 PetscInt row; 4519 row = i; 4520 if (rowindices) row = rowindices[i]; 4521 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4522 PetscInt col; 4523 col = Baij->j[count]; 4524 if (colindices) col = colindices[col]; 4525 v = Baij->a[count]; 4526 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4527 ++count; 4528 } 4529 } 4530 /* FIXME: set C's nonzerostate correctly. */ 4531 /* Assembly for C is necessary. */ 4532 C->preallocated = PETSC_TRUE; 4533 C->assembled = PETSC_TRUE; 4534 C->was_assembled = PETSC_FALSE; 4535 PetscFunctionReturn(0); 4536 } 4537 4538 4539 /* 4540 Special version for direct calls from Fortran 4541 */ 4542 #include <petsc/private/fortranimpl.h> 4543 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4544 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4545 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4546 #define matsetvaluesseqaij_ matsetvaluesseqaij 4547 #endif 4548 4549 /* Change these macros so can be used in void function */ 4550 #undef CHKERRQ 4551 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4552 #undef SETERRQ2 4553 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4554 #undef SETERRQ3 4555 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4556 4557 #undef __FUNCT__ 4558 #define __FUNCT__ "matsetvaluesseqaij_" 4559 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) 4560 { 4561 Mat A = *AA; 4562 PetscInt m = *mm, n = *nn; 4563 InsertMode is = *isis; 4564 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4565 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4566 PetscInt *imax,*ai,*ailen; 4567 PetscErrorCode ierr; 4568 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4569 MatScalar *ap,value,*aa; 4570 PetscBool ignorezeroentries = a->ignorezeroentries; 4571 PetscBool roworiented = a->roworiented; 4572 4573 PetscFunctionBegin; 4574 MatCheckPreallocated(A,1); 4575 imax = a->imax; 4576 ai = a->i; 4577 ailen = a->ilen; 4578 aj = a->j; 4579 aa = a->a; 4580 4581 for (k=0; k<m; k++) { /* loop over added rows */ 4582 row = im[k]; 4583 if (row < 0) continue; 4584 #if defined(PETSC_USE_DEBUG) 4585 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4586 #endif 4587 rp = aj + ai[row]; ap = aa + ai[row]; 4588 rmax = imax[row]; nrow = ailen[row]; 4589 low = 0; 4590 high = nrow; 4591 for (l=0; l<n; l++) { /* loop over added columns */ 4592 if (in[l] < 0) continue; 4593 #if defined(PETSC_USE_DEBUG) 4594 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4595 #endif 4596 col = in[l]; 4597 if (roworiented) value = v[l + k*n]; 4598 else value = v[k + l*m]; 4599 4600 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4601 4602 if (col <= lastcol) low = 0; 4603 else high = nrow; 4604 lastcol = col; 4605 while (high-low > 5) { 4606 t = (low+high)/2; 4607 if (rp[t] > col) high = t; 4608 else low = t; 4609 } 4610 for (i=low; i<high; i++) { 4611 if (rp[i] > col) break; 4612 if (rp[i] == col) { 4613 if (is == ADD_VALUES) ap[i] += value; 4614 else ap[i] = value; 4615 goto noinsert; 4616 } 4617 } 4618 if (value == 0.0 && ignorezeroentries) goto noinsert; 4619 if (nonew == 1) goto noinsert; 4620 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4621 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4622 N = nrow++ - 1; a->nz++; high++; 4623 /* shift up all the later entries in this row */ 4624 for (ii=N; ii>=i; ii--) { 4625 rp[ii+1] = rp[ii]; 4626 ap[ii+1] = ap[ii]; 4627 } 4628 rp[i] = col; 4629 ap[i] = value; 4630 A->nonzerostate++; 4631 noinsert:; 4632 low = i + 1; 4633 } 4634 ailen[row] = nrow; 4635 } 4636 PetscFunctionReturnVoid(); 4637 } 4638 4639