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