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