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