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