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