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