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