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