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]; /* omega in idiag */ 1704 } 1705 } 1706 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 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 /* lower */ 1714 n = diag[i] - a->i[i]; 1715 idx = a->j + a->i[i]; 1716 v = a->a + a->i[i]; 1717 sum = b[i]; 1718 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1719 t[i] = sum; /* save application of the lower-triangular part */ 1720 /* upper */ 1721 n = a->i[i+1] - diag[i] - 1; 1722 idx = a->j + diag[i] + 1; 1723 v = a->a + diag[i] + 1; 1724 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1725 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1726 } 1727 xb = t; 1728 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1729 } else xb = b; 1730 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1731 for (i=m-1; i>=0; i--) { 1732 sum = xb[i]; 1733 if (xb == b) { 1734 /* whole matrix (no checkpointing available) */ 1735 n = a->i[i+1] - a->i[i]; 1736 idx = a->j + a->i[i]; 1737 v = a->a + a->i[i]; 1738 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1739 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1740 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1741 n = a->i[i+1] - diag[i] - 1; 1742 idx = a->j + diag[i] + 1; 1743 v = a->a + diag[i] + 1; 1744 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1745 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1746 } 1747 } 1748 if (xb == b) { 1749 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1750 } else { 1751 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1752 } 1753 } 1754 } 1755 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1756 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1757 PetscFunctionReturn(0); 1758 } 1759 1760 1761 #undef __FUNCT__ 1762 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1763 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1764 { 1765 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1766 1767 PetscFunctionBegin; 1768 info->block_size = 1.0; 1769 info->nz_allocated = (double)a->maxnz; 1770 info->nz_used = (double)a->nz; 1771 info->nz_unneeded = (double)(a->maxnz - a->nz); 1772 info->assemblies = (double)A->num_ass; 1773 info->mallocs = (double)A->info.mallocs; 1774 info->memory = ((PetscObject)A)->mem; 1775 if (A->factortype) { 1776 info->fill_ratio_given = A->info.fill_ratio_given; 1777 info->fill_ratio_needed = A->info.fill_ratio_needed; 1778 info->factor_mallocs = A->info.factor_mallocs; 1779 } else { 1780 info->fill_ratio_given = 0; 1781 info->fill_ratio_needed = 0; 1782 info->factor_mallocs = 0; 1783 } 1784 PetscFunctionReturn(0); 1785 } 1786 1787 #undef __FUNCT__ 1788 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1789 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1790 { 1791 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1792 PetscInt i,m = A->rmap->n - 1,d = 0; 1793 PetscErrorCode ierr; 1794 const PetscScalar *xx; 1795 PetscScalar *bb; 1796 PetscBool missing; 1797 1798 PetscFunctionBegin; 1799 if (x && b) { 1800 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1801 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1802 for (i=0; i<N; i++) { 1803 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1804 bb[rows[i]] = diag*xx[rows[i]]; 1805 } 1806 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1807 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1808 } 1809 1810 if (a->keepnonzeropattern) { 1811 for (i=0; i<N; i++) { 1812 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1813 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1814 } 1815 if (diag != 0.0) { 1816 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1817 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1818 for (i=0; i<N; i++) { 1819 a->a[a->diag[rows[i]]] = diag; 1820 } 1821 } 1822 A->same_nonzero = PETSC_TRUE; 1823 } else { 1824 if (diag != 0.0) { 1825 for (i=0; i<N; i++) { 1826 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1827 if (a->ilen[rows[i]] > 0) { 1828 a->ilen[rows[i]] = 1; 1829 a->a[a->i[rows[i]]] = diag; 1830 a->j[a->i[rows[i]]] = rows[i]; 1831 } else { /* in case row was completely empty */ 1832 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1833 } 1834 } 1835 } else { 1836 for (i=0; i<N; i++) { 1837 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1838 a->ilen[rows[i]] = 0; 1839 } 1840 } 1841 A->same_nonzero = PETSC_FALSE; 1842 } 1843 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1844 PetscFunctionReturn(0); 1845 } 1846 1847 #undef __FUNCT__ 1848 #define __FUNCT__ "MatZeroRowsColumns_SeqAIJ" 1849 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1850 { 1851 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1852 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1853 PetscErrorCode ierr; 1854 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 1855 const PetscScalar *xx; 1856 PetscScalar *bb; 1857 1858 PetscFunctionBegin; 1859 if (x && b) { 1860 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1861 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1862 vecs = PETSC_TRUE; 1863 } 1864 ierr = PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);CHKERRQ(ierr); 1865 ierr = PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));CHKERRQ(ierr); 1866 for (i=0; i<N; i++) { 1867 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1868 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1869 1870 zeroed[rows[i]] = PETSC_TRUE; 1871 } 1872 for (i=0; i<A->rmap->n; i++) { 1873 if (!zeroed[i]) { 1874 for (j=a->i[i]; j<a->i[i+1]; j++) { 1875 if (zeroed[a->j[j]]) { 1876 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 1877 a->a[j] = 0.0; 1878 } 1879 } 1880 } else if (vecs) bb[i] = diag*xx[i]; 1881 } 1882 if (x && b) { 1883 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1884 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1885 } 1886 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1887 if (diag != 0.0) { 1888 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1889 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1890 for (i=0; i<N; i++) { 1891 a->a[a->diag[rows[i]]] = diag; 1892 } 1893 } 1894 A->same_nonzero = PETSC_TRUE; 1895 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1896 PetscFunctionReturn(0); 1897 } 1898 1899 #undef __FUNCT__ 1900 #define __FUNCT__ "MatGetRow_SeqAIJ" 1901 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1902 { 1903 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1904 PetscInt *itmp; 1905 1906 PetscFunctionBegin; 1907 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1908 1909 *nz = a->i[row+1] - a->i[row]; 1910 if (v) *v = a->a + a->i[row]; 1911 if (idx) { 1912 itmp = a->j + a->i[row]; 1913 if (*nz) *idx = itmp; 1914 else *idx = 0; 1915 } 1916 PetscFunctionReturn(0); 1917 } 1918 1919 /* remove this function? */ 1920 #undef __FUNCT__ 1921 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1922 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1923 { 1924 PetscFunctionBegin; 1925 PetscFunctionReturn(0); 1926 } 1927 1928 #undef __FUNCT__ 1929 #define __FUNCT__ "MatNorm_SeqAIJ" 1930 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1931 { 1932 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1933 MatScalar *v = a->a; 1934 PetscReal sum = 0.0; 1935 PetscErrorCode ierr; 1936 PetscInt i,j; 1937 1938 PetscFunctionBegin; 1939 if (type == NORM_FROBENIUS) { 1940 for (i=0; i<a->nz; i++) { 1941 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1942 } 1943 *nrm = PetscSqrtReal(sum); 1944 } else if (type == NORM_1) { 1945 PetscReal *tmp; 1946 PetscInt *jj = a->j; 1947 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1948 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1949 *nrm = 0.0; 1950 for (j=0; j<a->nz; j++) { 1951 tmp[*jj++] += PetscAbsScalar(*v); v++; 1952 } 1953 for (j=0; j<A->cmap->n; j++) { 1954 if (tmp[j] > *nrm) *nrm = tmp[j]; 1955 } 1956 ierr = PetscFree(tmp);CHKERRQ(ierr); 1957 } else if (type == NORM_INFINITY) { 1958 *nrm = 0.0; 1959 for (j=0; j<A->rmap->n; j++) { 1960 v = a->a + a->i[j]; 1961 sum = 0.0; 1962 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1963 sum += PetscAbsScalar(*v); v++; 1964 } 1965 if (sum > *nrm) *nrm = sum; 1966 } 1967 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1968 PetscFunctionReturn(0); 1969 } 1970 1971 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 1972 #undef __FUNCT__ 1973 #define __FUNCT__ "MatTransposeSymbolic_SeqAIJ" 1974 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 1975 { 1976 PetscErrorCode ierr; 1977 PetscInt i,j,anzj; 1978 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 1979 PetscInt an=A->cmap->N,am=A->rmap->N; 1980 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 1981 1982 PetscFunctionBegin; 1983 /* Allocate space for symbolic transpose info and work array */ 1984 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 1985 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 1986 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 1987 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 1988 1989 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 1990 /* Note: offset by 1 for fast conversion into csr format. */ 1991 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 1992 /* Form ati for csr format of A^T. */ 1993 for (i=0;i<an;i++) ati[i+1] += ati[i]; 1994 1995 /* Copy ati into atfill so we have locations of the next free space in atj */ 1996 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 1997 1998 /* Walk through A row-wise and mark nonzero entries of A^T. */ 1999 for (i=0;i<am;i++) { 2000 anzj = ai[i+1] - ai[i]; 2001 for (j=0;j<anzj;j++) { 2002 atj[atfill[*aj]] = i; 2003 atfill[*aj++] += 1; 2004 } 2005 } 2006 2007 /* Clean up temporary space and complete requests. */ 2008 ierr = PetscFree(atfill);CHKERRQ(ierr); 2009 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2010 2011 (*B)->rmap->bs = A->cmap->bs; 2012 (*B)->cmap->bs = A->rmap->bs; 2013 2014 b = (Mat_SeqAIJ*)((*B)->data); 2015 b->free_a = PETSC_FALSE; 2016 b->free_ij = PETSC_TRUE; 2017 b->nonew = 0; 2018 PetscFunctionReturn(0); 2019 } 2020 2021 #undef __FUNCT__ 2022 #define __FUNCT__ "MatTranspose_SeqAIJ" 2023 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2024 { 2025 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2026 Mat C; 2027 PetscErrorCode ierr; 2028 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2029 MatScalar *array = a->a; 2030 2031 PetscFunctionBegin; 2032 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"); 2033 2034 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 2035 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 2036 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 2037 2038 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2039 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2040 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2041 ierr = MatSetBlockSizes(C,A->cmap->bs,A->rmap->bs);CHKERRQ(ierr); 2042 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2043 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2044 ierr = PetscFree(col);CHKERRQ(ierr); 2045 } else { 2046 C = *B; 2047 } 2048 2049 for (i=0; i<m; i++) { 2050 len = ai[i+1]-ai[i]; 2051 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2052 array += len; 2053 aj += len; 2054 } 2055 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2056 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2057 2058 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 2059 *B = C; 2060 } else { 2061 ierr = MatHeaderMerge(A,C);CHKERRQ(ierr); 2062 } 2063 PetscFunctionReturn(0); 2064 } 2065 2066 #undef __FUNCT__ 2067 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 2068 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2069 { 2070 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2071 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2072 MatScalar *va,*vb; 2073 PetscErrorCode ierr; 2074 PetscInt ma,na,mb,nb, i; 2075 2076 PetscFunctionBegin; 2077 bij = (Mat_SeqAIJ*) B->data; 2078 2079 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2080 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2081 if (ma!=nb || na!=mb) { 2082 *f = PETSC_FALSE; 2083 PetscFunctionReturn(0); 2084 } 2085 aii = aij->i; bii = bij->i; 2086 adx = aij->j; bdx = bij->j; 2087 va = aij->a; vb = bij->a; 2088 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 2089 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 2090 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2091 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2092 2093 *f = PETSC_TRUE; 2094 for (i=0; i<ma; i++) { 2095 while (aptr[i]<aii[i+1]) { 2096 PetscInt idc,idr; 2097 PetscScalar vc,vr; 2098 /* column/row index/value */ 2099 idc = adx[aptr[i]]; 2100 idr = bdx[bptr[idc]]; 2101 vc = va[aptr[i]]; 2102 vr = vb[bptr[idc]]; 2103 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2104 *f = PETSC_FALSE; 2105 goto done; 2106 } else { 2107 aptr[i]++; 2108 if (B || i!=idc) bptr[idc]++; 2109 } 2110 } 2111 } 2112 done: 2113 ierr = PetscFree(aptr);CHKERRQ(ierr); 2114 ierr = PetscFree(bptr);CHKERRQ(ierr); 2115 PetscFunctionReturn(0); 2116 } 2117 2118 #undef __FUNCT__ 2119 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 2120 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2121 { 2122 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2123 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2124 MatScalar *va,*vb; 2125 PetscErrorCode ierr; 2126 PetscInt ma,na,mb,nb, i; 2127 2128 PetscFunctionBegin; 2129 bij = (Mat_SeqAIJ*) B->data; 2130 2131 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2132 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2133 if (ma!=nb || na!=mb) { 2134 *f = PETSC_FALSE; 2135 PetscFunctionReturn(0); 2136 } 2137 aii = aij->i; bii = bij->i; 2138 adx = aij->j; bdx = bij->j; 2139 va = aij->a; vb = bij->a; 2140 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 2141 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 2142 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2143 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2144 2145 *f = PETSC_TRUE; 2146 for (i=0; i<ma; i++) { 2147 while (aptr[i]<aii[i+1]) { 2148 PetscInt idc,idr; 2149 PetscScalar vc,vr; 2150 /* column/row index/value */ 2151 idc = adx[aptr[i]]; 2152 idr = bdx[bptr[idc]]; 2153 vc = va[aptr[i]]; 2154 vr = vb[bptr[idc]]; 2155 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2156 *f = PETSC_FALSE; 2157 goto done; 2158 } else { 2159 aptr[i]++; 2160 if (B || i!=idc) bptr[idc]++; 2161 } 2162 } 2163 } 2164 done: 2165 ierr = PetscFree(aptr);CHKERRQ(ierr); 2166 ierr = PetscFree(bptr);CHKERRQ(ierr); 2167 PetscFunctionReturn(0); 2168 } 2169 2170 #undef __FUNCT__ 2171 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 2172 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2173 { 2174 PetscErrorCode ierr; 2175 2176 PetscFunctionBegin; 2177 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2178 PetscFunctionReturn(0); 2179 } 2180 2181 #undef __FUNCT__ 2182 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 2183 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2184 { 2185 PetscErrorCode ierr; 2186 2187 PetscFunctionBegin; 2188 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2189 PetscFunctionReturn(0); 2190 } 2191 2192 #undef __FUNCT__ 2193 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 2194 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2195 { 2196 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2197 PetscScalar *l,*r,x; 2198 MatScalar *v; 2199 PetscErrorCode ierr; 2200 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 2201 2202 PetscFunctionBegin; 2203 if (ll) { 2204 /* The local size is used so that VecMPI can be passed to this routine 2205 by MatDiagonalScale_MPIAIJ */ 2206 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2207 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2208 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 2209 v = a->a; 2210 for (i=0; i<m; i++) { 2211 x = l[i]; 2212 M = a->i[i+1] - a->i[i]; 2213 for (j=0; j<M; j++) (*v++) *= x; 2214 } 2215 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 2216 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2217 } 2218 if (rr) { 2219 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2220 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2221 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 2222 v = a->a; jj = a->j; 2223 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2224 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 2225 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2226 } 2227 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2228 PetscFunctionReturn(0); 2229 } 2230 2231 #undef __FUNCT__ 2232 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 2233 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2234 { 2235 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2236 PetscErrorCode ierr; 2237 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2238 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2239 const PetscInt *irow,*icol; 2240 PetscInt nrows,ncols; 2241 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2242 MatScalar *a_new,*mat_a; 2243 Mat C; 2244 PetscBool stride,sorted; 2245 2246 PetscFunctionBegin; 2247 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 2248 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 2249 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 2250 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 2251 2252 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2253 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2254 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2255 2256 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2257 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2258 if (stride && step == 1) { 2259 /* special case of contiguous rows */ 2260 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 2261 /* loop over new rows determining lens and starting points */ 2262 for (i=0; i<nrows; i++) { 2263 kstart = ai[irow[i]]; 2264 kend = kstart + ailen[irow[i]]; 2265 for (k=kstart; k<kend; k++) { 2266 if (aj[k] >= first) { 2267 starts[i] = k; 2268 break; 2269 } 2270 } 2271 sum = 0; 2272 while (k < kend) { 2273 if (aj[k++] >= first+ncols) break; 2274 sum++; 2275 } 2276 lens[i] = sum; 2277 } 2278 /* create submatrix */ 2279 if (scall == MAT_REUSE_MATRIX) { 2280 PetscInt n_cols,n_rows; 2281 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2282 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2283 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2284 C = *B; 2285 } else { 2286 PetscInt rbs,cbs; 2287 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2288 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2289 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2290 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2291 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2292 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2293 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2294 } 2295 c = (Mat_SeqAIJ*)C->data; 2296 2297 /* loop over rows inserting into submatrix */ 2298 a_new = c->a; 2299 j_new = c->j; 2300 i_new = c->i; 2301 2302 for (i=0; i<nrows; i++) { 2303 ii = starts[i]; 2304 lensi = lens[i]; 2305 for (k=0; k<lensi; k++) { 2306 *j_new++ = aj[ii+k] - first; 2307 } 2308 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2309 a_new += lensi; 2310 i_new[i+1] = i_new[i] + lensi; 2311 c->ilen[i] = lensi; 2312 } 2313 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2314 } else { 2315 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2316 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 2317 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 2318 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2319 for (i=0; i<ncols; i++) { 2320 #if defined(PETSC_USE_DEBUG) 2321 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); 2322 #endif 2323 smap[icol[i]] = i+1; 2324 } 2325 2326 /* determine lens of each row */ 2327 for (i=0; i<nrows; i++) { 2328 kstart = ai[irow[i]]; 2329 kend = kstart + a->ilen[irow[i]]; 2330 lens[i] = 0; 2331 for (k=kstart; k<kend; k++) { 2332 if (smap[aj[k]]) { 2333 lens[i]++; 2334 } 2335 } 2336 } 2337 /* Create and fill new matrix */ 2338 if (scall == MAT_REUSE_MATRIX) { 2339 PetscBool equal; 2340 2341 c = (Mat_SeqAIJ*)((*B)->data); 2342 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2343 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2344 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2345 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2346 C = *B; 2347 } else { 2348 PetscInt rbs,cbs; 2349 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2350 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2351 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2352 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2353 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2354 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2355 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2356 } 2357 c = (Mat_SeqAIJ*)(C->data); 2358 for (i=0; i<nrows; i++) { 2359 row = irow[i]; 2360 kstart = ai[row]; 2361 kend = kstart + a->ilen[row]; 2362 mat_i = c->i[i]; 2363 mat_j = c->j + mat_i; 2364 mat_a = c->a + mat_i; 2365 mat_ilen = c->ilen + i; 2366 for (k=kstart; k<kend; k++) { 2367 if ((tcol=smap[a->j[k]])) { 2368 *mat_j++ = tcol - 1; 2369 *mat_a++ = a->a[k]; 2370 (*mat_ilen)++; 2371 2372 } 2373 } 2374 } 2375 /* Free work space */ 2376 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2377 ierr = PetscFree(smap);CHKERRQ(ierr); 2378 ierr = PetscFree(lens);CHKERRQ(ierr); 2379 } 2380 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2381 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2382 2383 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2384 *B = C; 2385 PetscFunctionReturn(0); 2386 } 2387 2388 #undef __FUNCT__ 2389 #define __FUNCT__ "MatGetMultiProcBlock_SeqAIJ" 2390 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2391 { 2392 PetscErrorCode ierr; 2393 Mat B; 2394 2395 PetscFunctionBegin; 2396 if (scall == MAT_INITIAL_MATRIX) { 2397 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2398 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2399 ierr = MatSetBlockSizes(B,mat->rmap->bs,mat->cmap->bs);CHKERRQ(ierr); 2400 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2401 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2402 *subMat = B; 2403 } else { 2404 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2405 } 2406 PetscFunctionReturn(0); 2407 } 2408 2409 #undef __FUNCT__ 2410 #define __FUNCT__ "MatILUFactor_SeqAIJ" 2411 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2412 { 2413 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2414 PetscErrorCode ierr; 2415 Mat outA; 2416 PetscBool row_identity,col_identity; 2417 2418 PetscFunctionBegin; 2419 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2420 2421 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2422 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2423 2424 outA = inA; 2425 outA->factortype = MAT_FACTOR_LU; 2426 2427 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2428 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2429 2430 a->row = row; 2431 2432 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2433 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2434 2435 a->col = col; 2436 2437 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2438 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2439 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2440 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2441 2442 if (!a->solve_work) { /* this matrix may have been factored before */ 2443 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 2444 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2445 } 2446 2447 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2448 if (row_identity && col_identity) { 2449 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2450 } else { 2451 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2452 } 2453 PetscFunctionReturn(0); 2454 } 2455 2456 #undef __FUNCT__ 2457 #define __FUNCT__ "MatScale_SeqAIJ" 2458 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2459 { 2460 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2461 PetscScalar oalpha = alpha; 2462 PetscErrorCode ierr; 2463 PetscBLASInt one = 1,bnz; 2464 2465 PetscFunctionBegin; 2466 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2467 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2468 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2469 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2470 PetscFunctionReturn(0); 2471 } 2472 2473 #undef __FUNCT__ 2474 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 2475 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2476 { 2477 PetscErrorCode ierr; 2478 PetscInt i; 2479 2480 PetscFunctionBegin; 2481 if (scall == MAT_INITIAL_MATRIX) { 2482 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 2483 } 2484 2485 for (i=0; i<n; i++) { 2486 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2487 } 2488 PetscFunctionReturn(0); 2489 } 2490 2491 #undef __FUNCT__ 2492 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 2493 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2494 { 2495 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2496 PetscErrorCode ierr; 2497 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2498 const PetscInt *idx; 2499 PetscInt start,end,*ai,*aj; 2500 PetscBT table; 2501 2502 PetscFunctionBegin; 2503 m = A->rmap->n; 2504 ai = a->i; 2505 aj = a->j; 2506 2507 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2508 2509 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 2510 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2511 2512 for (i=0; i<is_max; i++) { 2513 /* Initialize the two local arrays */ 2514 isz = 0; 2515 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2516 2517 /* Extract the indices, assume there can be duplicate entries */ 2518 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2519 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2520 2521 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2522 for (j=0; j<n; ++j) { 2523 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2524 } 2525 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2526 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2527 2528 k = 0; 2529 for (j=0; j<ov; j++) { /* for each overlap */ 2530 n = isz; 2531 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2532 row = nidx[k]; 2533 start = ai[row]; 2534 end = ai[row+1]; 2535 for (l = start; l<end; l++) { 2536 val = aj[l]; 2537 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2538 } 2539 } 2540 } 2541 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2542 } 2543 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2544 ierr = PetscFree(nidx);CHKERRQ(ierr); 2545 PetscFunctionReturn(0); 2546 } 2547 2548 /* -------------------------------------------------------------- */ 2549 #undef __FUNCT__ 2550 #define __FUNCT__ "MatPermute_SeqAIJ" 2551 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2552 { 2553 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2554 PetscErrorCode ierr; 2555 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2556 const PetscInt *row,*col; 2557 PetscInt *cnew,j,*lens; 2558 IS icolp,irowp; 2559 PetscInt *cwork = NULL; 2560 PetscScalar *vwork = NULL; 2561 2562 PetscFunctionBegin; 2563 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2564 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2565 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2566 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2567 2568 /* determine lengths of permuted rows */ 2569 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2570 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2571 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2572 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2573 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 2574 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2575 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2576 ierr = PetscFree(lens);CHKERRQ(ierr); 2577 2578 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2579 for (i=0; i<m; i++) { 2580 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2581 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2582 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2583 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2584 } 2585 ierr = PetscFree(cnew);CHKERRQ(ierr); 2586 2587 (*B)->assembled = PETSC_FALSE; 2588 2589 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2590 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2591 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2592 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2593 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2594 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2595 PetscFunctionReturn(0); 2596 } 2597 2598 #undef __FUNCT__ 2599 #define __FUNCT__ "MatCopy_SeqAIJ" 2600 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2601 { 2602 PetscErrorCode ierr; 2603 2604 PetscFunctionBegin; 2605 /* If the two matrices have the same copy implementation, use fast copy. */ 2606 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2607 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2608 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2609 2610 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"); 2611 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2612 } else { 2613 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2614 } 2615 PetscFunctionReturn(0); 2616 } 2617 2618 #undef __FUNCT__ 2619 #define __FUNCT__ "MatSetUp_SeqAIJ" 2620 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2621 { 2622 PetscErrorCode ierr; 2623 2624 PetscFunctionBegin; 2625 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2626 PetscFunctionReturn(0); 2627 } 2628 2629 #undef __FUNCT__ 2630 #define __FUNCT__ "MatSeqAIJGetArray_SeqAIJ" 2631 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2632 { 2633 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2634 2635 PetscFunctionBegin; 2636 *array = a->a; 2637 PetscFunctionReturn(0); 2638 } 2639 2640 #undef __FUNCT__ 2641 #define __FUNCT__ "MatSeqAIJRestoreArray_SeqAIJ" 2642 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2643 { 2644 PetscFunctionBegin; 2645 PetscFunctionReturn(0); 2646 } 2647 2648 #undef __FUNCT__ 2649 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2650 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2651 { 2652 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f; 2653 PetscErrorCode ierr; 2654 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow; 2655 PetscScalar dx,*y,*xx,*w3_array; 2656 PetscScalar *vscale_array; 2657 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2658 Vec w1,w2,w3; 2659 void *fctx = coloring->fctx; 2660 PetscBool flg = PETSC_FALSE; 2661 2662 PetscFunctionBegin; 2663 if (!coloring->w1) { 2664 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2665 ierr = PetscLogObjectParent((PetscObject)coloring,(PetscObject)coloring->w1);CHKERRQ(ierr); 2666 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2667 ierr = PetscLogObjectParent((PetscObject)coloring,(PetscObject)coloring->w2);CHKERRQ(ierr); 2668 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2669 ierr = PetscLogObjectParent((PetscObject)coloring,(PetscObject)coloring->w3);CHKERRQ(ierr); 2670 } 2671 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2672 2673 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2674 ierr = PetscOptionsGetBool(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,NULL);CHKERRQ(ierr); 2675 if (flg) { 2676 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2677 } else { 2678 PetscBool assembled; 2679 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2680 if (assembled) { 2681 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2682 } 2683 } 2684 2685 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2686 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2687 2688 if (!coloring->fset) { 2689 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2690 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2691 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2692 } else { 2693 coloring->fset = PETSC_FALSE; 2694 } 2695 2696 /* 2697 Compute all the scale factors and share with other processors 2698 */ 2699 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2700 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2701 for (k=0; k<coloring->ncolors; k++) { 2702 /* 2703 Loop over each column associated with color adding the 2704 perturbation to the vector w3. 2705 */ 2706 for (l=0; l<coloring->ncolumns[k]; l++) { 2707 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2708 dx = xx[col]; 2709 if (dx == 0.0) dx = 1.0; 2710 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2711 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2712 dx *= epsilon; 2713 vscale_array[col] = 1.0/dx; 2714 } 2715 } 2716 vscale_array = vscale_array + start; 2717 2718 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2719 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2720 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2721 2722 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2723 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2724 2725 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2726 else vscaleforrow = coloring->columnsforrow; 2727 2728 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2729 /* 2730 Loop over each color 2731 */ 2732 for (k=0; k<coloring->ncolors; k++) { 2733 coloring->currentcolor = k; 2734 2735 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2736 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2737 /* 2738 Loop over each column associated with color adding the 2739 perturbation to the vector w3. 2740 */ 2741 for (l=0; l<coloring->ncolumns[k]; l++) { 2742 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2743 dx = xx[col]; 2744 if (dx == 0.0) dx = 1.0; 2745 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2746 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2747 dx *= epsilon; 2748 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2749 w3_array[col] += dx; 2750 } 2751 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2752 2753 /* 2754 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2755 */ 2756 2757 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2758 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2759 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2760 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2761 2762 /* 2763 Loop over rows of vector, putting results into Jacobian matrix 2764 */ 2765 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2766 for (l=0; l<coloring->nrows[k]; l++) { 2767 row = coloring->rows[k][l]; 2768 col = coloring->columnsforrow[k][l]; 2769 y[row] *= vscale_array[vscaleforrow[k][l]]; 2770 srow = row + start; 2771 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2772 } 2773 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2774 } 2775 coloring->currentcolor = k; 2776 2777 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2778 xx = xx + start; 2779 ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2780 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2781 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2782 PetscFunctionReturn(0); 2783 } 2784 2785 /* 2786 Computes the number of nonzeros per row needed for preallocation when X and Y 2787 have different nonzero structure. 2788 */ 2789 #undef __FUNCT__ 2790 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2791 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2792 { 2793 PetscInt i,m=Y->rmap->N; 2794 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2795 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2796 const PetscInt *xi = x->i,*yi = y->i; 2797 2798 PetscFunctionBegin; 2799 /* Set the number of nonzeros in the new matrix */ 2800 for (i=0; i<m; i++) { 2801 PetscInt j,k,nzx = xi[i+1] - xi[i],nzy = yi[i+1] - yi[i]; 2802 const PetscInt *xj = x->j+xi[i],*yj = y->j+yi[i]; 2803 nnz[i] = 0; 2804 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2805 for (; k<nzy && yj[k]<xj[j]; k++) nnz[i]++; /* Catch up to X */ 2806 if (k<nzy && yj[k]==xj[j]) k++; /* Skip duplicate */ 2807 nnz[i]++; 2808 } 2809 for (; k<nzy; k++) nnz[i]++; 2810 } 2811 PetscFunctionReturn(0); 2812 } 2813 2814 #undef __FUNCT__ 2815 #define __FUNCT__ "MatAXPY_SeqAIJ" 2816 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2817 { 2818 PetscErrorCode ierr; 2819 PetscInt i; 2820 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2821 PetscBLASInt one=1,bnz; 2822 2823 PetscFunctionBegin; 2824 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2825 if (str == SAME_NONZERO_PATTERN) { 2826 PetscScalar alpha = a; 2827 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2828 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2829 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2830 if (y->xtoy && y->XtoY != X) { 2831 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2832 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 2833 } 2834 if (!y->xtoy) { /* get xtoy */ 2835 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,NULL, y->i,y->j,NULL, &y->xtoy);CHKERRQ(ierr); 2836 y->XtoY = X; 2837 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2838 } 2839 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2840 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); 2841 } else { 2842 Mat B; 2843 PetscInt *nnz; 2844 ierr = PetscMalloc(Y->rmap->N*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 2845 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2846 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2847 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2848 ierr = MatSetBlockSizes(B,Y->rmap->bs,Y->cmap->bs);CHKERRQ(ierr); 2849 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2850 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2851 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2852 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2853 ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr); 2854 ierr = PetscFree(nnz);CHKERRQ(ierr); 2855 } 2856 PetscFunctionReturn(0); 2857 } 2858 2859 #undef __FUNCT__ 2860 #define __FUNCT__ "MatConjugate_SeqAIJ" 2861 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2862 { 2863 #if defined(PETSC_USE_COMPLEX) 2864 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2865 PetscInt i,nz; 2866 PetscScalar *a; 2867 2868 PetscFunctionBegin; 2869 nz = aij->nz; 2870 a = aij->a; 2871 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2872 #else 2873 PetscFunctionBegin; 2874 #endif 2875 PetscFunctionReturn(0); 2876 } 2877 2878 #undef __FUNCT__ 2879 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2880 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2881 { 2882 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2883 PetscErrorCode ierr; 2884 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2885 PetscReal atmp; 2886 PetscScalar *x; 2887 MatScalar *aa; 2888 2889 PetscFunctionBegin; 2890 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2891 aa = a->a; 2892 ai = a->i; 2893 aj = a->j; 2894 2895 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2896 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2897 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2898 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2899 for (i=0; i<m; i++) { 2900 ncols = ai[1] - ai[0]; ai++; 2901 x[i] = 0.0; 2902 for (j=0; j<ncols; j++) { 2903 atmp = PetscAbsScalar(*aa); 2904 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2905 aa++; aj++; 2906 } 2907 } 2908 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2909 PetscFunctionReturn(0); 2910 } 2911 2912 #undef __FUNCT__ 2913 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2914 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2915 { 2916 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2917 PetscErrorCode ierr; 2918 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2919 PetscScalar *x; 2920 MatScalar *aa; 2921 2922 PetscFunctionBegin; 2923 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2924 aa = a->a; 2925 ai = a->i; 2926 aj = a->j; 2927 2928 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2929 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2930 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2931 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2932 for (i=0; i<m; i++) { 2933 ncols = ai[1] - ai[0]; ai++; 2934 if (ncols == A->cmap->n) { /* row is dense */ 2935 x[i] = *aa; if (idx) idx[i] = 0; 2936 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2937 x[i] = 0.0; 2938 if (idx) { 2939 idx[i] = 0; /* in case ncols is zero */ 2940 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2941 if (aj[j] > j) { 2942 idx[i] = j; 2943 break; 2944 } 2945 } 2946 } 2947 } 2948 for (j=0; j<ncols; j++) { 2949 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2950 aa++; aj++; 2951 } 2952 } 2953 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2954 PetscFunctionReturn(0); 2955 } 2956 2957 #undef __FUNCT__ 2958 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2959 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2960 { 2961 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2962 PetscErrorCode ierr; 2963 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2964 PetscReal atmp; 2965 PetscScalar *x; 2966 MatScalar *aa; 2967 2968 PetscFunctionBegin; 2969 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2970 aa = a->a; 2971 ai = a->i; 2972 aj = a->j; 2973 2974 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2975 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2976 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2977 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); 2978 for (i=0; i<m; i++) { 2979 ncols = ai[1] - ai[0]; ai++; 2980 if (ncols) { 2981 /* Get first nonzero */ 2982 for (j = 0; j < ncols; j++) { 2983 atmp = PetscAbsScalar(aa[j]); 2984 if (atmp > 1.0e-12) { 2985 x[i] = atmp; 2986 if (idx) idx[i] = aj[j]; 2987 break; 2988 } 2989 } 2990 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2991 } else { 2992 x[i] = 0.0; if (idx) idx[i] = 0; 2993 } 2994 for (j = 0; j < ncols; j++) { 2995 atmp = PetscAbsScalar(*aa); 2996 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2997 aa++; aj++; 2998 } 2999 } 3000 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3001 PetscFunctionReturn(0); 3002 } 3003 3004 #undef __FUNCT__ 3005 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 3006 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3007 { 3008 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3009 PetscErrorCode ierr; 3010 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3011 PetscScalar *x; 3012 MatScalar *aa; 3013 3014 PetscFunctionBegin; 3015 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3016 aa = a->a; 3017 ai = a->i; 3018 aj = a->j; 3019 3020 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3021 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3022 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3023 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3024 for (i=0; i<m; i++) { 3025 ncols = ai[1] - ai[0]; ai++; 3026 if (ncols == A->cmap->n) { /* row is dense */ 3027 x[i] = *aa; if (idx) idx[i] = 0; 3028 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3029 x[i] = 0.0; 3030 if (idx) { /* find first implicit 0.0 in the row */ 3031 idx[i] = 0; /* in case ncols is zero */ 3032 for (j=0; j<ncols; j++) { 3033 if (aj[j] > j) { 3034 idx[i] = j; 3035 break; 3036 } 3037 } 3038 } 3039 } 3040 for (j=0; j<ncols; j++) { 3041 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3042 aa++; aj++; 3043 } 3044 } 3045 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3046 PetscFunctionReturn(0); 3047 } 3048 3049 #include <petscblaslapack.h> 3050 #include <petsc-private/kernels/blockinvert.h> 3051 3052 #undef __FUNCT__ 3053 #define __FUNCT__ "MatInvertBlockDiagonal_SeqAIJ" 3054 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3055 { 3056 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3057 PetscErrorCode ierr; 3058 PetscInt i,bs = A->rmap->bs,mbs = A->rmap->n/A->rmap->bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3059 MatScalar *diag,work[25],*v_work; 3060 PetscReal shift = 0.0; 3061 3062 PetscFunctionBegin; 3063 if (a->ibdiagvalid) { 3064 if (values) *values = a->ibdiag; 3065 PetscFunctionReturn(0); 3066 } 3067 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 3068 if (!a->ibdiag) { 3069 ierr = PetscMalloc(bs2*mbs*sizeof(PetscScalar),&a->ibdiag);CHKERRQ(ierr); 3070 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 3071 } 3072 diag = a->ibdiag; 3073 if (values) *values = a->ibdiag; 3074 /* factor and invert each block */ 3075 switch (bs) { 3076 case 1: 3077 for (i=0; i<mbs; i++) { 3078 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3079 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3080 } 3081 break; 3082 case 2: 3083 for (i=0; i<mbs; i++) { 3084 ij[0] = 2*i; ij[1] = 2*i + 1; 3085 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3086 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift);CHKERRQ(ierr); 3087 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3088 diag += 4; 3089 } 3090 break; 3091 case 3: 3092 for (i=0; i<mbs; i++) { 3093 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3094 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3095 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift);CHKERRQ(ierr); 3096 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3097 diag += 9; 3098 } 3099 break; 3100 case 4: 3101 for (i=0; i<mbs; i++) { 3102 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3103 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3104 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift);CHKERRQ(ierr); 3105 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3106 diag += 16; 3107 } 3108 break; 3109 case 5: 3110 for (i=0; i<mbs; i++) { 3111 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3112 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3113 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift);CHKERRQ(ierr); 3114 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3115 diag += 25; 3116 } 3117 break; 3118 case 6: 3119 for (i=0; i<mbs; i++) { 3120 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; 3121 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3122 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift);CHKERRQ(ierr); 3123 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3124 diag += 36; 3125 } 3126 break; 3127 case 7: 3128 for (i=0; i<mbs; i++) { 3129 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; 3130 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3131 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift);CHKERRQ(ierr); 3132 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3133 diag += 49; 3134 } 3135 break; 3136 default: 3137 ierr = PetscMalloc3(bs,MatScalar,&v_work,bs,PetscInt,&v_pivots,bs,PetscInt,&IJ);CHKERRQ(ierr); 3138 for (i=0; i<mbs; i++) { 3139 for (j=0; j<bs; j++) { 3140 IJ[j] = bs*i + j; 3141 } 3142 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3143 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);CHKERRQ(ierr); 3144 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3145 diag += bs2; 3146 } 3147 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3148 } 3149 a->ibdiagvalid = PETSC_TRUE; 3150 PetscFunctionReturn(0); 3151 } 3152 3153 #undef __FUNCT__ 3154 #define __FUNCT__ "MatSetRandom_SeqAIJ" 3155 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3156 { 3157 PetscErrorCode ierr; 3158 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3159 PetscScalar a; 3160 PetscInt m,n,i,j,col; 3161 3162 PetscFunctionBegin; 3163 if (!x->assembled) { 3164 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3165 for (i=0; i<m; i++) { 3166 for (j=0; j<aij->imax[i]; j++) { 3167 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3168 col = (PetscInt)(n*PetscRealPart(a)); 3169 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3170 } 3171 } 3172 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3173 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3174 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3175 PetscFunctionReturn(0); 3176 } 3177 3178 extern PetscErrorCode MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 3179 /* -------------------------------------------------------------------*/ 3180 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3181 MatGetRow_SeqAIJ, 3182 MatRestoreRow_SeqAIJ, 3183 MatMult_SeqAIJ, 3184 /* 4*/ MatMultAdd_SeqAIJ, 3185 MatMultTranspose_SeqAIJ, 3186 MatMultTransposeAdd_SeqAIJ, 3187 0, 3188 0, 3189 0, 3190 /* 10*/ 0, 3191 MatLUFactor_SeqAIJ, 3192 0, 3193 MatSOR_SeqAIJ, 3194 MatTranspose_SeqAIJ, 3195 /*1 5*/ MatGetInfo_SeqAIJ, 3196 MatEqual_SeqAIJ, 3197 MatGetDiagonal_SeqAIJ, 3198 MatDiagonalScale_SeqAIJ, 3199 MatNorm_SeqAIJ, 3200 /* 20*/ 0, 3201 MatAssemblyEnd_SeqAIJ, 3202 MatSetOption_SeqAIJ, 3203 MatZeroEntries_SeqAIJ, 3204 /* 24*/ MatZeroRows_SeqAIJ, 3205 0, 3206 0, 3207 0, 3208 0, 3209 /* 29*/ MatSetUp_SeqAIJ, 3210 0, 3211 0, 3212 0, 3213 0, 3214 /* 34*/ MatDuplicate_SeqAIJ, 3215 0, 3216 0, 3217 MatILUFactor_SeqAIJ, 3218 0, 3219 /* 39*/ MatAXPY_SeqAIJ, 3220 MatGetSubMatrices_SeqAIJ, 3221 MatIncreaseOverlap_SeqAIJ, 3222 MatGetValues_SeqAIJ, 3223 MatCopy_SeqAIJ, 3224 /* 44*/ MatGetRowMax_SeqAIJ, 3225 MatScale_SeqAIJ, 3226 0, 3227 MatDiagonalSet_SeqAIJ, 3228 MatZeroRowsColumns_SeqAIJ, 3229 /* 49*/ MatSetRandom_SeqAIJ, 3230 MatGetRowIJ_SeqAIJ, 3231 MatRestoreRowIJ_SeqAIJ, 3232 MatGetColumnIJ_SeqAIJ, 3233 MatRestoreColumnIJ_SeqAIJ, 3234 /* 54*/ MatFDColoringCreate_SeqAIJ, 3235 0, 3236 0, 3237 MatPermute_SeqAIJ, 3238 0, 3239 /* 59*/ 0, 3240 MatDestroy_SeqAIJ, 3241 MatView_SeqAIJ, 3242 0, 3243 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3244 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3245 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3246 0, 3247 0, 3248 0, 3249 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3250 MatGetRowMinAbs_SeqAIJ, 3251 0, 3252 MatSetColoring_SeqAIJ, 3253 0, 3254 /* 74*/ MatSetValuesAdifor_SeqAIJ, 3255 MatFDColoringApply_AIJ, 3256 0, 3257 0, 3258 0, 3259 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3260 0, 3261 0, 3262 0, 3263 MatLoad_SeqAIJ, 3264 /* 84*/ MatIsSymmetric_SeqAIJ, 3265 MatIsHermitian_SeqAIJ, 3266 0, 3267 0, 3268 0, 3269 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3270 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3271 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3272 MatPtAP_SeqAIJ_SeqAIJ, 3273 MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy, 3274 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3275 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3276 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3277 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3278 0, 3279 /* 99*/ 0, 3280 0, 3281 0, 3282 MatConjugate_SeqAIJ, 3283 0, 3284 /*104*/ MatSetValuesRow_SeqAIJ, 3285 MatRealPart_SeqAIJ, 3286 MatImaginaryPart_SeqAIJ, 3287 0, 3288 0, 3289 /*109*/ MatMatSolve_SeqAIJ, 3290 0, 3291 MatGetRowMin_SeqAIJ, 3292 0, 3293 MatMissingDiagonal_SeqAIJ, 3294 /*114*/ 0, 3295 0, 3296 0, 3297 0, 3298 0, 3299 /*119*/ 0, 3300 0, 3301 0, 3302 0, 3303 MatGetMultiProcBlock_SeqAIJ, 3304 /*124*/ MatFindNonzeroRows_SeqAIJ, 3305 MatGetColumnNorms_SeqAIJ, 3306 MatInvertBlockDiagonal_SeqAIJ, 3307 0, 3308 0, 3309 /*129*/ 0, 3310 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3311 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3312 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3313 MatTransposeColoringCreate_SeqAIJ, 3314 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3315 MatTransColoringApplyDenToSp_SeqAIJ, 3316 MatRARt_SeqAIJ_SeqAIJ, 3317 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3318 MatRARtNumeric_SeqAIJ_SeqAIJ, 3319 /*139*/0, 3320 0 3321 }; 3322 3323 #undef __FUNCT__ 3324 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3325 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3326 { 3327 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3328 PetscInt i,nz,n; 3329 3330 PetscFunctionBegin; 3331 nz = aij->maxnz; 3332 n = mat->rmap->n; 3333 for (i=0; i<nz; i++) { 3334 aij->j[i] = indices[i]; 3335 } 3336 aij->nz = nz; 3337 for (i=0; i<n; i++) { 3338 aij->ilen[i] = aij->imax[i]; 3339 } 3340 PetscFunctionReturn(0); 3341 } 3342 3343 #undef __FUNCT__ 3344 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3345 /*@ 3346 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3347 in the matrix. 3348 3349 Input Parameters: 3350 + mat - the SeqAIJ matrix 3351 - indices - the column indices 3352 3353 Level: advanced 3354 3355 Notes: 3356 This can be called if you have precomputed the nonzero structure of the 3357 matrix and want to provide it to the matrix object to improve the performance 3358 of the MatSetValues() operation. 3359 3360 You MUST have set the correct numbers of nonzeros per row in the call to 3361 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3362 3363 MUST be called before any calls to MatSetValues(); 3364 3365 The indices should start with zero, not one. 3366 3367 @*/ 3368 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3369 { 3370 PetscErrorCode ierr; 3371 3372 PetscFunctionBegin; 3373 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3374 PetscValidPointer(indices,2); 3375 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3376 PetscFunctionReturn(0); 3377 } 3378 3379 /* ----------------------------------------------------------------------------------------*/ 3380 3381 #undef __FUNCT__ 3382 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3383 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3384 { 3385 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3386 PetscErrorCode ierr; 3387 size_t nz = aij->i[mat->rmap->n]; 3388 3389 PetscFunctionBegin; 3390 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3391 3392 /* allocate space for values if not already there */ 3393 if (!aij->saved_values) { 3394 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 3395 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3396 } 3397 3398 /* copy values over */ 3399 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3400 PetscFunctionReturn(0); 3401 } 3402 3403 #undef __FUNCT__ 3404 #define __FUNCT__ "MatStoreValues" 3405 /*@ 3406 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3407 example, reuse of the linear part of a Jacobian, while recomputing the 3408 nonlinear portion. 3409 3410 Collect on Mat 3411 3412 Input Parameters: 3413 . mat - the matrix (currently only AIJ matrices support this option) 3414 3415 Level: advanced 3416 3417 Common Usage, with SNESSolve(): 3418 $ Create Jacobian matrix 3419 $ Set linear terms into matrix 3420 $ Apply boundary conditions to matrix, at this time matrix must have 3421 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3422 $ boundary conditions again will not change the nonzero structure 3423 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3424 $ ierr = MatStoreValues(mat); 3425 $ Call SNESSetJacobian() with matrix 3426 $ In your Jacobian routine 3427 $ ierr = MatRetrieveValues(mat); 3428 $ Set nonlinear terms in matrix 3429 3430 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3431 $ // build linear portion of Jacobian 3432 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3433 $ ierr = MatStoreValues(mat); 3434 $ loop over nonlinear iterations 3435 $ ierr = MatRetrieveValues(mat); 3436 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3437 $ // call MatAssemblyBegin/End() on matrix 3438 $ Solve linear system with Jacobian 3439 $ endloop 3440 3441 Notes: 3442 Matrix must already be assemblied before calling this routine 3443 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3444 calling this routine. 3445 3446 When this is called multiple times it overwrites the previous set of stored values 3447 and does not allocated additional space. 3448 3449 .seealso: MatRetrieveValues() 3450 3451 @*/ 3452 PetscErrorCode MatStoreValues(Mat mat) 3453 { 3454 PetscErrorCode ierr; 3455 3456 PetscFunctionBegin; 3457 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3458 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3459 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3460 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3461 PetscFunctionReturn(0); 3462 } 3463 3464 #undef __FUNCT__ 3465 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3466 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3467 { 3468 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3469 PetscErrorCode ierr; 3470 PetscInt nz = aij->i[mat->rmap->n]; 3471 3472 PetscFunctionBegin; 3473 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3474 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3475 /* copy values over */ 3476 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3477 PetscFunctionReturn(0); 3478 } 3479 3480 #undef __FUNCT__ 3481 #define __FUNCT__ "MatRetrieveValues" 3482 /*@ 3483 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3484 example, reuse of the linear part of a Jacobian, while recomputing the 3485 nonlinear portion. 3486 3487 Collect on Mat 3488 3489 Input Parameters: 3490 . mat - the matrix (currently on AIJ matrices support this option) 3491 3492 Level: advanced 3493 3494 .seealso: MatStoreValues() 3495 3496 @*/ 3497 PetscErrorCode MatRetrieveValues(Mat mat) 3498 { 3499 PetscErrorCode ierr; 3500 3501 PetscFunctionBegin; 3502 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3503 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3504 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3505 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3506 PetscFunctionReturn(0); 3507 } 3508 3509 3510 /* --------------------------------------------------------------------------------*/ 3511 #undef __FUNCT__ 3512 #define __FUNCT__ "MatCreateSeqAIJ" 3513 /*@C 3514 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3515 (the default parallel PETSc format). For good matrix assembly performance 3516 the user should preallocate the matrix storage by setting the parameter nz 3517 (or the array nnz). By setting these parameters accurately, performance 3518 during matrix assembly can be increased by more than a factor of 50. 3519 3520 Collective on MPI_Comm 3521 3522 Input Parameters: 3523 + comm - MPI communicator, set to PETSC_COMM_SELF 3524 . m - number of rows 3525 . n - number of columns 3526 . nz - number of nonzeros per row (same for all rows) 3527 - nnz - array containing the number of nonzeros in the various rows 3528 (possibly different for each row) or NULL 3529 3530 Output Parameter: 3531 . A - the matrix 3532 3533 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3534 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3535 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3536 3537 Notes: 3538 If nnz is given then nz is ignored 3539 3540 The AIJ format (also called the Yale sparse matrix format or 3541 compressed row storage), is fully compatible with standard Fortran 77 3542 storage. That is, the stored row and column indices can begin at 3543 either one (as in Fortran) or zero. See the users' manual for details. 3544 3545 Specify the preallocated storage with either nz or nnz (not both). 3546 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3547 allocation. For large problems you MUST preallocate memory or you 3548 will get TERRIBLE performance, see the users' manual chapter on matrices. 3549 3550 By default, this format uses inodes (identical nodes) when possible, to 3551 improve numerical efficiency of matrix-vector products and solves. We 3552 search for consecutive rows with the same nonzero structure, thereby 3553 reusing matrix information to achieve increased efficiency. 3554 3555 Options Database Keys: 3556 + -mat_no_inode - Do not use inodes 3557 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3558 3559 Level: intermediate 3560 3561 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3562 3563 @*/ 3564 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3565 { 3566 PetscErrorCode ierr; 3567 3568 PetscFunctionBegin; 3569 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3570 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3571 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3572 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3573 PetscFunctionReturn(0); 3574 } 3575 3576 #undef __FUNCT__ 3577 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3578 /*@C 3579 MatSeqAIJSetPreallocation - For good matrix assembly performance 3580 the user should preallocate the matrix storage by setting the parameter nz 3581 (or the array nnz). By setting these parameters accurately, performance 3582 during matrix assembly can be increased by more than a factor of 50. 3583 3584 Collective on MPI_Comm 3585 3586 Input Parameters: 3587 + B - The matrix-free 3588 . nz - number of nonzeros per row (same for all rows) 3589 - nnz - array containing the number of nonzeros in the various rows 3590 (possibly different for each row) or NULL 3591 3592 Notes: 3593 If nnz is given then nz is ignored 3594 3595 The AIJ format (also called the Yale sparse matrix format or 3596 compressed row storage), is fully compatible with standard Fortran 77 3597 storage. That is, the stored row and column indices can begin at 3598 either one (as in Fortran) or zero. See the users' manual for details. 3599 3600 Specify the preallocated storage with either nz or nnz (not both). 3601 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3602 allocation. For large problems you MUST preallocate memory or you 3603 will get TERRIBLE performance, see the users' manual chapter on matrices. 3604 3605 You can call MatGetInfo() to get information on how effective the preallocation was; 3606 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3607 You can also run with the option -info and look for messages with the string 3608 malloc in them to see if additional memory allocation was needed. 3609 3610 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3611 entries or columns indices 3612 3613 By default, this format uses inodes (identical nodes) when possible, to 3614 improve numerical efficiency of matrix-vector products and solves. We 3615 search for consecutive rows with the same nonzero structure, thereby 3616 reusing matrix information to achieve increased efficiency. 3617 3618 Options Database Keys: 3619 + -mat_no_inode - Do not use inodes 3620 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3621 - -mat_aij_oneindex - Internally use indexing starting at 1 3622 rather than 0. Note that when calling MatSetValues(), 3623 the user still MUST index entries starting at 0! 3624 3625 Level: intermediate 3626 3627 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3628 3629 @*/ 3630 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3631 { 3632 PetscErrorCode ierr; 3633 3634 PetscFunctionBegin; 3635 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3636 PetscValidType(B,1); 3637 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3638 PetscFunctionReturn(0); 3639 } 3640 3641 #undef __FUNCT__ 3642 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3643 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3644 { 3645 Mat_SeqAIJ *b; 3646 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3647 PetscErrorCode ierr; 3648 PetscInt i; 3649 3650 PetscFunctionBegin; 3651 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3652 if (nz == MAT_SKIP_ALLOCATION) { 3653 skipallocation = PETSC_TRUE; 3654 nz = 0; 3655 } 3656 3657 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3658 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3659 3660 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3661 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3662 if (nnz) { 3663 for (i=0; i<B->rmap->n; i++) { 3664 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]); 3665 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); 3666 } 3667 } 3668 3669 B->preallocated = PETSC_TRUE; 3670 3671 b = (Mat_SeqAIJ*)B->data; 3672 3673 if (!skipallocation) { 3674 if (!b->imax) { 3675 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3676 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3677 } 3678 if (!nnz) { 3679 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3680 else if (nz < 0) nz = 1; 3681 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3682 nz = nz*B->rmap->n; 3683 } else { 3684 nz = 0; 3685 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3686 } 3687 /* b->ilen will count nonzeros in each row so far. */ 3688 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3689 3690 /* allocate the matrix space */ 3691 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3692 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3693 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3694 b->i[0] = 0; 3695 for (i=1; i<B->rmap->n+1; i++) { 3696 b->i[i] = b->i[i-1] + b->imax[i-1]; 3697 } 3698 b->singlemalloc = PETSC_TRUE; 3699 b->free_a = PETSC_TRUE; 3700 b->free_ij = PETSC_TRUE; 3701 #if defined(PETSC_THREADCOMM_ACTIVE) 3702 ierr = MatZeroEntries_SeqAIJ(B);CHKERRQ(ierr); 3703 #endif 3704 } else { 3705 b->free_a = PETSC_FALSE; 3706 b->free_ij = PETSC_FALSE; 3707 } 3708 3709 b->nz = 0; 3710 b->maxnz = nz; 3711 B->info.nz_unneeded = (double)b->maxnz; 3712 if (realalloc) { 3713 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3714 } 3715 PetscFunctionReturn(0); 3716 } 3717 3718 #undef __FUNCT__ 3719 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3720 /*@ 3721 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3722 3723 Input Parameters: 3724 + B - the matrix 3725 . i - the indices into j for the start of each row (starts with zero) 3726 . j - the column indices for each row (starts with zero) these must be sorted for each row 3727 - v - optional values in the matrix 3728 3729 Level: developer 3730 3731 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3732 3733 .keywords: matrix, aij, compressed row, sparse, sequential 3734 3735 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3736 @*/ 3737 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3738 { 3739 PetscErrorCode ierr; 3740 3741 PetscFunctionBegin; 3742 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3743 PetscValidType(B,1); 3744 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3745 PetscFunctionReturn(0); 3746 } 3747 3748 #undef __FUNCT__ 3749 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3750 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3751 { 3752 PetscInt i; 3753 PetscInt m,n; 3754 PetscInt nz; 3755 PetscInt *nnz, nz_max = 0; 3756 PetscScalar *values; 3757 PetscErrorCode ierr; 3758 3759 PetscFunctionBegin; 3760 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3761 3762 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3763 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3764 3765 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3766 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3767 for (i = 0; i < m; i++) { 3768 nz = Ii[i+1]- Ii[i]; 3769 nz_max = PetscMax(nz_max, nz); 3770 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3771 nnz[i] = nz; 3772 } 3773 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3774 ierr = PetscFree(nnz);CHKERRQ(ierr); 3775 3776 if (v) { 3777 values = (PetscScalar*) v; 3778 } else { 3779 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3780 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3781 } 3782 3783 for (i = 0; i < m; i++) { 3784 nz = Ii[i+1] - Ii[i]; 3785 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3786 } 3787 3788 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3789 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3790 3791 if (!v) { 3792 ierr = PetscFree(values);CHKERRQ(ierr); 3793 } 3794 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3795 PetscFunctionReturn(0); 3796 } 3797 3798 #include <../src/mat/impls/dense/seq/dense.h> 3799 #include <petsc-private/kernels/petscaxpy.h> 3800 3801 #undef __FUNCT__ 3802 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3803 /* 3804 Computes (B'*A')' since computing B*A directly is untenable 3805 3806 n p p 3807 ( ) ( ) ( ) 3808 m ( A ) * n ( B ) = m ( C ) 3809 ( ) ( ) ( ) 3810 3811 */ 3812 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3813 { 3814 PetscErrorCode ierr; 3815 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3816 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3817 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3818 PetscInt i,n,m,q,p; 3819 const PetscInt *ii,*idx; 3820 const PetscScalar *b,*a,*a_q; 3821 PetscScalar *c,*c_q; 3822 3823 PetscFunctionBegin; 3824 m = A->rmap->n; 3825 n = A->cmap->n; 3826 p = B->cmap->n; 3827 a = sub_a->v; 3828 b = sub_b->a; 3829 c = sub_c->v; 3830 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3831 3832 ii = sub_b->i; 3833 idx = sub_b->j; 3834 for (i=0; i<n; i++) { 3835 q = ii[i+1] - ii[i]; 3836 while (q-->0) { 3837 c_q = c + m*(*idx); 3838 a_q = a + m*i; 3839 PetscKernelAXPY(c_q,*b,a_q,m); 3840 idx++; 3841 b++; 3842 } 3843 } 3844 PetscFunctionReturn(0); 3845 } 3846 3847 #undef __FUNCT__ 3848 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3849 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3850 { 3851 PetscErrorCode ierr; 3852 PetscInt m=A->rmap->n,n=B->cmap->n; 3853 Mat Cmat; 3854 3855 PetscFunctionBegin; 3856 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); 3857 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3858 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3859 ierr = MatSetBlockSizes(Cmat,A->rmap->bs,B->cmap->bs);CHKERRQ(ierr); 3860 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3861 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3862 3863 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3864 3865 *C = Cmat; 3866 PetscFunctionReturn(0); 3867 } 3868 3869 /* ----------------------------------------------------------------*/ 3870 #undef __FUNCT__ 3871 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3872 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3873 { 3874 PetscErrorCode ierr; 3875 3876 PetscFunctionBegin; 3877 if (scall == MAT_INITIAL_MATRIX) { 3878 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3879 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3880 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3881 } 3882 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3883 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3884 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3885 PetscFunctionReturn(0); 3886 } 3887 3888 3889 /*MC 3890 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3891 based on compressed sparse row format. 3892 3893 Options Database Keys: 3894 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3895 3896 Level: beginner 3897 3898 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3899 M*/ 3900 3901 /*MC 3902 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3903 3904 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3905 and MATMPIAIJ otherwise. As a result, for single process communicators, 3906 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3907 for communicators controlling multiple processes. It is recommended that you call both of 3908 the above preallocation routines for simplicity. 3909 3910 Options Database Keys: 3911 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3912 3913 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3914 enough exist. 3915 3916 Level: beginner 3917 3918 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3919 M*/ 3920 3921 /*MC 3922 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3923 3924 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3925 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3926 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3927 for communicators controlling multiple processes. It is recommended that you call both of 3928 the above preallocation routines for simplicity. 3929 3930 Options Database Keys: 3931 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3932 3933 Level: beginner 3934 3935 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3936 M*/ 3937 3938 #if defined(PETSC_HAVE_PASTIX) 3939 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3940 #endif 3941 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3942 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat*); 3943 #endif 3944 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3945 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3946 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3947 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool*); 3948 #if defined(PETSC_HAVE_MUMPS) 3949 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3950 #endif 3951 #if defined(PETSC_HAVE_SUPERLU) 3952 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3953 #endif 3954 #if defined(PETSC_HAVE_SUPERLU_DIST) 3955 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3956 #endif 3957 #if defined(PETSC_HAVE_UMFPACK) 3958 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3959 #endif 3960 #if defined(PETSC_HAVE_CHOLMOD) 3961 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3962 #endif 3963 #if defined(PETSC_HAVE_LUSOL) 3964 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3965 #endif 3966 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3967 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3968 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3969 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3970 #endif 3971 #if defined(PETSC_HAVE_CLIQUE) 3972 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_clique(Mat,MatFactorType,Mat*); 3973 #endif 3974 3975 3976 #undef __FUNCT__ 3977 #define __FUNCT__ "MatSeqAIJGetArray" 3978 /*@C 3979 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 3980 3981 Not Collective 3982 3983 Input Parameter: 3984 . mat - a MATSEQDENSE matrix 3985 3986 Output Parameter: 3987 . array - pointer to the data 3988 3989 Level: intermediate 3990 3991 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3992 @*/ 3993 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3994 { 3995 PetscErrorCode ierr; 3996 3997 PetscFunctionBegin; 3998 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3999 PetscFunctionReturn(0); 4000 } 4001 4002 #undef __FUNCT__ 4003 #define __FUNCT__ "MatSeqAIJRestoreArray" 4004 /*@C 4005 MatSeqAIJRestoreArray - returns access to the array where the data for a SeqSeqAIJ matrix is stored obtained by MatSeqAIJGetArray() 4006 4007 Not Collective 4008 4009 Input Parameters: 4010 . mat - a MATSEQDENSE matrix 4011 . array - pointer to the data 4012 4013 Level: intermediate 4014 4015 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4016 @*/ 4017 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4018 { 4019 PetscErrorCode ierr; 4020 4021 PetscFunctionBegin; 4022 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4023 PetscFunctionReturn(0); 4024 } 4025 4026 #undef __FUNCT__ 4027 #define __FUNCT__ "MatCreate_SeqAIJ" 4028 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4029 { 4030 Mat_SeqAIJ *b; 4031 PetscErrorCode ierr; 4032 PetscMPIInt size; 4033 4034 PetscFunctionBegin; 4035 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4036 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4037 4038 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 4039 4040 B->data = (void*)b; 4041 4042 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4043 4044 b->row = 0; 4045 b->col = 0; 4046 b->icol = 0; 4047 b->reallocs = 0; 4048 b->ignorezeroentries = PETSC_FALSE; 4049 b->roworiented = PETSC_TRUE; 4050 b->nonew = 0; 4051 b->diag = 0; 4052 b->solve_work = 0; 4053 B->spptr = 0; 4054 b->saved_values = 0; 4055 b->idiag = 0; 4056 b->mdiag = 0; 4057 b->ssor_work = 0; 4058 b->omega = 1.0; 4059 b->fshift = 0.0; 4060 b->idiagvalid = PETSC_FALSE; 4061 b->ibdiagvalid = PETSC_FALSE; 4062 b->keepnonzeropattern = PETSC_FALSE; 4063 b->xtoy = 0; 4064 b->XtoY = 0; 4065 B->same_nonzero = PETSC_FALSE; 4066 4067 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4068 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4069 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4070 4071 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4072 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_matlab_C",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 4073 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4074 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4075 #endif 4076 #if defined(PETSC_HAVE_PASTIX) 4077 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_pastix_C",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 4078 #endif 4079 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 4080 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_essl_C",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 4081 #endif 4082 #if defined(PETSC_HAVE_SUPERLU) 4083 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_C",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 4084 #endif 4085 #if defined(PETSC_HAVE_SUPERLU_DIST) 4086 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_dist_C",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 4087 #endif 4088 #if defined(PETSC_HAVE_MUMPS) 4089 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_mumps_C",MatGetFactor_aij_mumps);CHKERRQ(ierr); 4090 #endif 4091 #if defined(PETSC_HAVE_UMFPACK) 4092 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_umfpack_C",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 4093 #endif 4094 #if defined(PETSC_HAVE_CHOLMOD) 4095 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_cholmod_C",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 4096 #endif 4097 #if defined(PETSC_HAVE_LUSOL) 4098 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_lusol_C",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 4099 #endif 4100 #if defined(PETSC_HAVE_CLIQUE) 4101 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_clique_C",MatGetFactor_aij_clique);CHKERRQ(ierr); 4102 #endif 4103 4104 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 4105 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactorAvailable_petsc_C",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 4106 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_bas_C",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 4107 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4108 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4109 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4110 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4111 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4112 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4113 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4114 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4115 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4116 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4117 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4118 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4119 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4120 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4121 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4122 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4123 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4124 PetscFunctionReturn(0); 4125 } 4126 4127 #undef __FUNCT__ 4128 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 4129 /* 4130 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4131 */ 4132 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4133 { 4134 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4135 PetscErrorCode ierr; 4136 PetscInt i,m = A->rmap->n; 4137 4138 PetscFunctionBegin; 4139 c = (Mat_SeqAIJ*)C->data; 4140 4141 C->factortype = A->factortype; 4142 c->row = 0; 4143 c->col = 0; 4144 c->icol = 0; 4145 c->reallocs = 0; 4146 4147 C->assembled = PETSC_TRUE; 4148 4149 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4150 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4151 4152 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 4153 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4154 for (i=0; i<m; i++) { 4155 c->imax[i] = a->imax[i]; 4156 c->ilen[i] = a->ilen[i]; 4157 } 4158 4159 /* allocate the matrix space */ 4160 if (mallocmatspace) { 4161 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 4162 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4163 4164 c->singlemalloc = PETSC_TRUE; 4165 4166 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4167 if (m > 0) { 4168 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4169 if (cpvalues == MAT_COPY_VALUES) { 4170 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4171 } else { 4172 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4173 } 4174 } 4175 } 4176 4177 c->ignorezeroentries = a->ignorezeroentries; 4178 c->roworiented = a->roworiented; 4179 c->nonew = a->nonew; 4180 if (a->diag) { 4181 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 4182 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4183 for (i=0; i<m; i++) { 4184 c->diag[i] = a->diag[i]; 4185 } 4186 } else c->diag = 0; 4187 4188 c->solve_work = 0; 4189 c->saved_values = 0; 4190 c->idiag = 0; 4191 c->ssor_work = 0; 4192 c->keepnonzeropattern = a->keepnonzeropattern; 4193 c->free_a = PETSC_TRUE; 4194 c->free_ij = PETSC_TRUE; 4195 c->xtoy = 0; 4196 c->XtoY = 0; 4197 4198 c->rmax = a->rmax; 4199 c->nz = a->nz; 4200 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4201 C->preallocated = PETSC_TRUE; 4202 4203 c->compressedrow.use = a->compressedrow.use; 4204 c->compressedrow.nrows = a->compressedrow.nrows; 4205 c->compressedrow.check = a->compressedrow.check; 4206 if (a->compressedrow.use) { 4207 i = a->compressedrow.nrows; 4208 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 4209 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4210 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4211 } else { 4212 c->compressedrow.use = PETSC_FALSE; 4213 c->compressedrow.i = NULL; 4214 c->compressedrow.rindex = NULL; 4215 } 4216 C->same_nonzero = A->same_nonzero; 4217 4218 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4219 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4220 PetscFunctionReturn(0); 4221 } 4222 4223 #undef __FUNCT__ 4224 #define __FUNCT__ "MatDuplicate_SeqAIJ" 4225 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4226 { 4227 PetscErrorCode ierr; 4228 4229 PetscFunctionBegin; 4230 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4231 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4232 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 4233 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4234 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4235 PetscFunctionReturn(0); 4236 } 4237 4238 #undef __FUNCT__ 4239 #define __FUNCT__ "MatLoad_SeqAIJ" 4240 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4241 { 4242 Mat_SeqAIJ *a; 4243 PetscErrorCode ierr; 4244 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4245 int fd; 4246 PetscMPIInt size; 4247 MPI_Comm comm; 4248 PetscInt bs = 1; 4249 4250 PetscFunctionBegin; 4251 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4252 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4253 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4254 4255 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4256 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4257 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4258 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 4259 4260 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4261 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4262 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4263 M = header[1]; N = header[2]; nz = header[3]; 4264 4265 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4266 4267 /* read in row lengths */ 4268 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 4269 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4270 4271 /* check if sum of rowlengths is same as nz */ 4272 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4273 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); 4274 4275 /* set global size if not set already*/ 4276 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4277 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4278 } else { 4279 /* if sizes and type are already set, check if the vector global sizes are correct */ 4280 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4281 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4282 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4283 } 4284 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); 4285 } 4286 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4287 a = (Mat_SeqAIJ*)newMat->data; 4288 4289 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4290 4291 /* read in nonzero values */ 4292 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4293 4294 /* set matrix "i" values */ 4295 a->i[0] = 0; 4296 for (i=1; i<= M; i++) { 4297 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4298 a->ilen[i-1] = rowlengths[i-1]; 4299 } 4300 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4301 4302 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4303 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4304 PetscFunctionReturn(0); 4305 } 4306 4307 #undef __FUNCT__ 4308 #define __FUNCT__ "MatEqual_SeqAIJ" 4309 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4310 { 4311 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4312 PetscErrorCode ierr; 4313 #if defined(PETSC_USE_COMPLEX) 4314 PetscInt k; 4315 #endif 4316 4317 PetscFunctionBegin; 4318 /* If the matrix dimensions are not equal,or no of nonzeros */ 4319 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4320 *flg = PETSC_FALSE; 4321 PetscFunctionReturn(0); 4322 } 4323 4324 /* if the a->i are the same */ 4325 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4326 if (!*flg) PetscFunctionReturn(0); 4327 4328 /* if a->j are the same */ 4329 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4330 if (!*flg) PetscFunctionReturn(0); 4331 4332 /* if a->a are the same */ 4333 #if defined(PETSC_USE_COMPLEX) 4334 for (k=0; k<a->nz; k++) { 4335 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4336 *flg = PETSC_FALSE; 4337 PetscFunctionReturn(0); 4338 } 4339 } 4340 #else 4341 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4342 #endif 4343 PetscFunctionReturn(0); 4344 } 4345 4346 #undef __FUNCT__ 4347 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4348 /*@ 4349 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4350 provided by the user. 4351 4352 Collective on MPI_Comm 4353 4354 Input Parameters: 4355 + comm - must be an MPI communicator of size 1 4356 . m - number of rows 4357 . n - number of columns 4358 . i - row indices 4359 . j - column indices 4360 - a - matrix values 4361 4362 Output Parameter: 4363 . mat - the matrix 4364 4365 Level: intermediate 4366 4367 Notes: 4368 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4369 once the matrix is destroyed and not before 4370 4371 You cannot set new nonzero locations into this matrix, that will generate an error. 4372 4373 The i and j indices are 0 based 4374 4375 The format which is used for the sparse matrix input, is equivalent to a 4376 row-major ordering.. i.e for the following matrix, the input data expected is 4377 as shown: 4378 4379 1 0 0 4380 2 0 3 4381 4 5 6 4382 4383 i = {0,1,3,6} [size = nrow+1 = 3+1] 4384 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4385 v = {1,2,3,4,5,6} [size = nz = 6] 4386 4387 4388 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4389 4390 @*/ 4391 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 4392 { 4393 PetscErrorCode ierr; 4394 PetscInt ii; 4395 Mat_SeqAIJ *aij; 4396 #if defined(PETSC_USE_DEBUG) 4397 PetscInt jj; 4398 #endif 4399 4400 PetscFunctionBegin; 4401 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4402 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4403 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4404 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4405 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4406 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4407 aij = (Mat_SeqAIJ*)(*mat)->data; 4408 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 4409 4410 aij->i = i; 4411 aij->j = j; 4412 aij->a = a; 4413 aij->singlemalloc = PETSC_FALSE; 4414 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4415 aij->free_a = PETSC_FALSE; 4416 aij->free_ij = PETSC_FALSE; 4417 4418 for (ii=0; ii<m; ii++) { 4419 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4420 #if defined(PETSC_USE_DEBUG) 4421 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]); 4422 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4423 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); 4424 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); 4425 } 4426 #endif 4427 } 4428 #if defined(PETSC_USE_DEBUG) 4429 for (ii=0; ii<aij->i[m]; ii++) { 4430 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 4431 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]); 4432 } 4433 #endif 4434 4435 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4436 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4437 PetscFunctionReturn(0); 4438 } 4439 #undef __FUNCT__ 4440 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4441 /*@C 4442 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4443 provided by the user. 4444 4445 Collective on MPI_Comm 4446 4447 Input Parameters: 4448 + comm - must be an MPI communicator of size 1 4449 . m - number of rows 4450 . n - number of columns 4451 . i - row indices 4452 . j - column indices 4453 . a - matrix values 4454 . nz - number of nonzeros 4455 - idx - 0 or 1 based 4456 4457 Output Parameter: 4458 . mat - the matrix 4459 4460 Level: intermediate 4461 4462 Notes: 4463 The i and j indices are 0 based 4464 4465 The format which is used for the sparse matrix input, is equivalent to a 4466 row-major ordering.. i.e for the following matrix, the input data expected is 4467 as shown: 4468 4469 1 0 0 4470 2 0 3 4471 4 5 6 4472 4473 i = {0,1,1,2,2,2} 4474 j = {0,0,2,0,1,2} 4475 v = {1,2,3,4,5,6} 4476 4477 4478 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4479 4480 @*/ 4481 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4482 { 4483 PetscErrorCode ierr; 4484 PetscInt ii, *nnz, one = 1,row,col; 4485 4486 4487 PetscFunctionBegin; 4488 ierr = PetscMalloc(m*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 4489 ierr = PetscMemzero(nnz,m*sizeof(PetscInt));CHKERRQ(ierr); 4490 for (ii = 0; ii < nz; ii++) { 4491 nnz[i[ii]] += 1; 4492 } 4493 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4494 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4495 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4496 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4497 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4498 for (ii = 0; ii < nz; ii++) { 4499 if (idx) { 4500 row = i[ii] - 1; 4501 col = j[ii] - 1; 4502 } else { 4503 row = i[ii]; 4504 col = j[ii]; 4505 } 4506 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4507 } 4508 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4509 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4510 ierr = PetscFree(nnz);CHKERRQ(ierr); 4511 PetscFunctionReturn(0); 4512 } 4513 4514 #undef __FUNCT__ 4515 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4516 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4517 { 4518 PetscErrorCode ierr; 4519 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4520 4521 PetscFunctionBegin; 4522 if (coloring->ctype == IS_COLORING_GLOBAL) { 4523 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4524 a->coloring = coloring; 4525 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4526 PetscInt i,*larray; 4527 ISColoring ocoloring; 4528 ISColoringValue *colors; 4529 4530 /* set coloring for diagonal portion */ 4531 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 4532 for (i=0; i<A->cmap->n; i++) larray[i] = i; 4533 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,NULL,larray);CHKERRQ(ierr); 4534 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 4535 for (i=0; i<A->cmap->n; i++) colors[i] = coloring->colors[larray[i]]; 4536 ierr = PetscFree(larray);CHKERRQ(ierr); 4537 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4538 a->coloring = ocoloring; 4539 } 4540 PetscFunctionReturn(0); 4541 } 4542 4543 #undef __FUNCT__ 4544 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4545 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4546 { 4547 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4548 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4549 MatScalar *v = a->a; 4550 PetscScalar *values = (PetscScalar*)advalues; 4551 ISColoringValue *color; 4552 4553 PetscFunctionBegin; 4554 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4555 color = a->coloring->colors; 4556 /* loop over rows */ 4557 for (i=0; i<m; i++) { 4558 nz = ii[i+1] - ii[i]; 4559 /* loop over columns putting computed value into matrix */ 4560 for (j=0; j<nz; j++) *v++ = values[color[*jj++]]; 4561 values += nl; /* jump to next row of derivatives */ 4562 } 4563 PetscFunctionReturn(0); 4564 } 4565 4566 #undef __FUNCT__ 4567 #define __FUNCT__ "MatSeqAIJInvalidateDiagonal" 4568 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4569 { 4570 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4571 PetscErrorCode ierr; 4572 4573 PetscFunctionBegin; 4574 a->idiagvalid = PETSC_FALSE; 4575 a->ibdiagvalid = PETSC_FALSE; 4576 4577 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4578 PetscFunctionReturn(0); 4579 } 4580 4581 /* 4582 Special version for direct calls from Fortran 4583 */ 4584 #include <petsc-private/fortranimpl.h> 4585 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4586 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4587 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4588 #define matsetvaluesseqaij_ matsetvaluesseqaij 4589 #endif 4590 4591 /* Change these macros so can be used in void function */ 4592 #undef CHKERRQ 4593 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4594 #undef SETERRQ2 4595 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4596 #undef SETERRQ3 4597 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4598 4599 #undef __FUNCT__ 4600 #define __FUNCT__ "matsetvaluesseqaij_" 4601 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) 4602 { 4603 Mat A = *AA; 4604 PetscInt m = *mm, n = *nn; 4605 InsertMode is = *isis; 4606 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4607 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4608 PetscInt *imax,*ai,*ailen; 4609 PetscErrorCode ierr; 4610 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4611 MatScalar *ap,value,*aa; 4612 PetscBool ignorezeroentries = a->ignorezeroentries; 4613 PetscBool roworiented = a->roworiented; 4614 4615 PetscFunctionBegin; 4616 MatCheckPreallocated(A,1); 4617 imax = a->imax; 4618 ai = a->i; 4619 ailen = a->ilen; 4620 aj = a->j; 4621 aa = a->a; 4622 4623 for (k=0; k<m; k++) { /* loop over added rows */ 4624 row = im[k]; 4625 if (row < 0) continue; 4626 #if defined(PETSC_USE_DEBUG) 4627 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4628 #endif 4629 rp = aj + ai[row]; ap = aa + ai[row]; 4630 rmax = imax[row]; nrow = ailen[row]; 4631 low = 0; 4632 high = nrow; 4633 for (l=0; l<n; l++) { /* loop over added columns */ 4634 if (in[l] < 0) continue; 4635 #if defined(PETSC_USE_DEBUG) 4636 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4637 #endif 4638 col = in[l]; 4639 if (roworiented) value = v[l + k*n]; 4640 else value = v[k + l*m]; 4641 4642 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4643 4644 if (col <= lastcol) low = 0; 4645 else high = nrow; 4646 lastcol = col; 4647 while (high-low > 5) { 4648 t = (low+high)/2; 4649 if (rp[t] > col) high = t; 4650 else low = t; 4651 } 4652 for (i=low; i<high; i++) { 4653 if (rp[i] > col) break; 4654 if (rp[i] == col) { 4655 if (is == ADD_VALUES) ap[i] += value; 4656 else ap[i] = value; 4657 goto noinsert; 4658 } 4659 } 4660 if (value == 0.0 && ignorezeroentries) goto noinsert; 4661 if (nonew == 1) goto noinsert; 4662 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4663 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4664 N = nrow++ - 1; a->nz++; high++; 4665 /* shift up all the later entries in this row */ 4666 for (ii=N; ii>=i; ii--) { 4667 rp[ii+1] = rp[ii]; 4668 ap[ii+1] = ap[ii]; 4669 } 4670 rp[i] = col; 4671 ap[i] = value; 4672 noinsert:; 4673 low = i + 1; 4674 } 4675 ailen[row] = nrow; 4676 } 4677 A->same_nonzero = PETSC_FALSE; 4678 PetscFunctionReturnVoid(); 4679 } 4680 4681 4682