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