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