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