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