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