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