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