SALSA Analysis Modules
|
00001 /*! \file simple.c \ingroup categories 00002 \brief Norm-like properties 00003 00004 \section simple Simple (normlike) quantities 00005 00006 The "simple" module contains simple, normlike quantities, 00007 In general, 00008 these are properties that can be calculated in time proportional to the 00009 number of nonzeros of the matrix. Note that this excludes the 2-norm. 00010 00011 \section Usage 00012 00013 Activate this module with 00014 00015 PetscErrorCode RegisterSimpleModules(); 00016 00017 Compute these elements with 00018 00019 ComputeQuantity("simple",<element>,A,(AnalysisItem*)&res,&flg); 00020 00021 Available elements are: 00022 - "trace" : Trace, sum of diagonal elements; see Trace(). 00023 - "trace-abs" : Trace Absolute, sum of absolute values of diagonal elements; 00024 see TraceAbs(). 00025 - "norm1" : 1-Norm, maximum column sum of absolute element sizes; 00026 see norm1(). 00027 - "normInf" : Infinity-Norm, maximum row sum of absolute element sizes; 00028 see normInf(). 00029 - "normF" : Frobenius-norm, sqrt of sum of elements squared; 00030 see normF(). 00031 - "diagonal-dominance" : least positive or most negative value of 00032 diagonal element minus sum of absolute off-diagonal elements; 00033 see diagonaldominance(). 00034 - "symmetry-snorm" : infinity norm of symmetric part; see SymmetrySNorm(). 00035 This element and all following are sequential; see \ref options; 00036 - "symmetry-anorm" : infinity norm of anti-symmetric part; see SymmetryANorm(). 00037 - "symmetry-fsnorm" : Frobenius norm of symmetric part; see 00038 SymmetryFSNorm(). 00039 - "symmetry-fanorm" : Frobenius norm of anti-symmetric part; 00040 see SymmetryFANorm(). 00041 - "n-struct-unsymm" : number of structurally unsymmetric elements; 00042 see NUnstruct(). 00043 (this is calculated here, but declared as an element of the \ref structure 00044 category.) 00045 */ 00046 00047 #include <stdlib.h> 00048 #include "anamod.h" 00049 #include "anamodsalsamodules.h" 00050 #include "petscerror.h" 00051 #include "petscmat.h" 00052 #include "src/mat/impls/aij/seq/aij.h" 00053 00054 #undef __FUNCT__ 00055 #define __FUNCT__ "computetrace" 00056 /*! 00057 This is the computational routine for trace calculation. 00058 It computes both the trace and the absolute trace, that is, using 00059 absolute values of matrix elements. 00060 00061 A Petsc vector is created and destroyed. 00062 */ 00063 static PetscErrorCode computetrace(Mat A) 00064 { 00065 MPI_Comm comm; Vec D; PetscScalar *d; 00066 PetscTruth has; PetscScalar vsum,gvsum; PetscReal asum,gasum; 00067 int local_size,i,id; 00068 PetscErrorCode ierr; 00069 00070 PetscFunctionBegin; 00071 00072 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00073 ierr = MatGetLocalSize(A,&local_size,PETSC_NULL); CHKERRQ(ierr); 00074 ierr = VecCreateMPI(comm,local_size,PETSC_DECIDE,&D); CHKERRQ(ierr); 00075 ierr = MatGetDiagonal(A,D); CHKERRQ(ierr); 00076 ierr = VecGetArray(D,&d); CHKERRQ(ierr); 00077 vsum = 0; asum = 0; 00078 for (i=0; i<local_size; i++) { 00079 vsum += d[i]; asum += PetscAbsScalar(d[i]); 00080 } 00081 ierr = VecRestoreArray(D,&d); CHKERRQ(ierr); 00082 ierr = VecDestroy(D); CHKERRQ(ierr); 00083 ierr = MPI_Allreduce 00084 ((void*)&vsum,(void*)&gvsum,1,MPIU_SCALAR,MPI_SUM,comm); CHKERRQ(ierr); 00085 ierr = MPI_Allreduce 00086 ((void*)&asum,(void*)&gasum,1,MPI_DOUBLE,MPI_SUM,comm); CHKERRQ(ierr); 00087 00088 ierr = GetDataID("simple","trace",&id,&has); CHKERRQ(ierr); 00089 HASTOEXIST(has); 00090 ierr = PetscObjectComposedDataSetScalar 00091 ((PetscObject)A,id,gvsum); CHKERRQ(ierr); 00092 ierr = GetDataID("simple","trace-abs",&id,&has); CHKERRQ(ierr); 00093 HASTOEXIST(has); 00094 ierr = PetscObjectComposedDataSetReal 00095 ((PetscObject)A,id,gasum); CHKERRQ(ierr); 00096 00097 PetscFunctionReturn(0); 00098 } 00099 00100 #undef __FUNCT__ 00101 #define __FUNCT__ "Trace" 00102 /*! Compute and store the trace of a matrix; 00103 the computation is done in computetrace(). 00104 */ 00105 static PetscErrorCode Trace 00106 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00107 { 00108 Mat A = (Mat)prob; 00109 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00110 PetscFunctionBegin; 00111 ierr = GetDataID("simple","trace",&id,&has); CHKERRQ(ierr); 00112 HASTOEXIST(has); 00113 ierr = PetscObjectComposedDataGetScalar 00114 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00115 if (*flg) { 00116 rv->r = v; 00117 } else { 00118 ierr = computetrace(A); CHKERRQ(ierr); 00119 ierr = PetscObjectComposedDataGetScalar 00120 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00121 if (*flg) { 00122 rv->r = v; 00123 } else SETERRQ(1,"Should have computed trace by now"); 00124 } 00125 00126 PetscFunctionReturn(0); 00127 } 00128 00129 #undef __FUNCT__ 00130 #define __FUNCT__ "TraceAbs" 00131 /*! Compute and store the absolute trace of a matrix; 00132 the computation is done in computetrace(). 00133 */ 00134 static PetscErrorCode TraceAbs 00135 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00136 { 00137 Mat A = (Mat)prob; 00138 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00139 PetscFunctionBegin; 00140 ierr = GetDataID("simple","trace-abs",&id,&has); CHKERRQ(ierr); 00141 HASTOEXIST(has); 00142 ierr = PetscObjectComposedDataGetReal 00143 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00144 if (*flg) rv->r = v; 00145 else { 00146 ierr = computetrace(A); CHKERRQ(ierr); 00147 ierr = PetscObjectComposedDataGetReal 00148 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00149 if (*flg) rv->r = v; 00150 } 00151 00152 PetscFunctionReturn(0); 00153 } 00154 00155 #undef __FUNCT__ 00156 #define __FUNCT__ "norm1" 00157 /*! Compute the 1-norm of a matrix. */ 00158 static PetscErrorCode norm1 00159 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00160 { 00161 Mat A = (Mat)prob; 00162 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00163 PetscFunctionBegin; 00164 00165 ierr = GetDataID("simple","norm1",&id,&has); CHKERRQ(ierr); 00166 HASTOEXIST(has); 00167 ierr = PetscObjectComposedDataGetReal 00168 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00169 if (*flg) rv->r = v; 00170 else { 00171 ierr = MatNorm(A,NORM_1,&v); CHKERRQ(ierr); 00172 ierr = PetscObjectComposedDataSetReal 00173 ((PetscObject)A,id,v); CHKERRQ(ierr); 00174 *flg = PETSC_TRUE; 00175 rv->r = v; 00176 } 00177 00178 *flg = PETSC_TRUE; 00179 PetscFunctionReturn(0); 00180 } 00181 00182 #undef __FUNCT__ 00183 #define __FUNCT__ "normF" 00184 /*! Compute the Frobenius norm of a matrix. */ 00185 static PetscErrorCode normF 00186 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00187 { 00188 Mat A = (Mat)prob; 00189 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00190 PetscFunctionBegin; 00191 00192 ierr = GetDataID("simple","normF",&id,&has); CHKERRQ(ierr); 00193 HASTOEXIST(has); 00194 ierr = PetscObjectComposedDataGetReal 00195 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00196 if (*flg) rv->r = v; 00197 else { 00198 ierr = MatNorm(A,NORM_FROBENIUS,&v); CHKERRQ(ierr); 00199 ierr = PetscObjectComposedDataSetReal 00200 ((PetscObject)A,id,v); CHKERRQ(ierr); 00201 *flg = PETSC_TRUE; 00202 rv->r = v; 00203 } 00204 00205 *flg = PETSC_TRUE; 00206 PetscFunctionReturn(0); 00207 } 00208 00209 #undef __FUNCT__ 00210 #define __FUNCT__ "normInf" 00211 /*! Compute the infinity-norm of a matrix. */ 00212 static PetscErrorCode normInf 00213 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00214 { 00215 Mat A = (Mat)prob; 00216 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00217 PetscFunctionBegin; 00218 00219 ierr = GetDataID("simple","normInf",&id,&has); CHKERRQ(ierr); 00220 HASTOEXIST(has); 00221 ierr = PetscObjectComposedDataGetReal 00222 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00223 if (*flg) rv->r = v; 00224 else { 00225 ierr = MatNorm(A,NORM_INFINITY,&v); CHKERRQ(ierr); 00226 ierr = PetscObjectComposedDataSetReal 00227 ((PetscObject)A,id,v); CHKERRQ(ierr); 00228 *flg = PETSC_TRUE; 00229 rv->r = v; 00230 } 00231 00232 *flg = PETSC_TRUE; 00233 PetscFunctionReturn(0); 00234 } 00235 00236 #undef __FUNCT__ 00237 #define __FUNCT__ "compute_dd" 00238 /*! 00239 Compute the diagonal dominance of a matrix: 00240 minimum value of diagonal element (not absolute!) minus off-diagonal 00241 elements absolute. 00242 \latexonly 00243 \[ \min_i a_{ii}-\sum_{j\not=i}|a_{ij}| \] 00244 \endlatexonly 00245 */ 00246 static PetscErrorCode compute_dd(Mat A,PetscTruth *flg) 00247 { 00248 PetscReal r,gr; 00249 MPI_Comm comm; 00250 int first,last,m,row; PetscErrorCode ierr; 00251 00252 PetscFunctionBegin; 00253 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00254 ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr); 00255 m = last-first; 00256 for (row=first; row<last; row++) { 00257 int col,ncols; const int *cols; 00258 const PetscScalar *vals; PetscReal rr=0.; 00259 ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr); 00260 for (col=0; col<ncols; col++) { 00261 PetscReal v = vals[col], vr = PetscRealPart(v), va = PetscAbsScalar(v); 00262 if (cols[col]==row) rr +=vr; else rr -= va; 00263 } 00264 if (row==first) r = rr; else r = PetscMin(r,rr); 00265 ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr); 00266 } 00267 MPI_Allreduce((void*)&r,(void*)&gr,1,MPIU_REAL,MPI_MIN,comm); 00268 { 00269 int id; PetscTruth has; 00270 ierr = GetDataID("simple","diagonal-dominance",&id,&has); CHKERRQ(ierr); 00271 HASTOEXIST(has); 00272 ierr = PetscObjectComposedDataSetReal 00273 ((PetscObject)A,id,gr); CHKERRQ(ierr); 00274 } 00275 *flg = PETSC_TRUE; 00276 PetscFunctionReturn(0); 00277 } 00278 00279 static PetscErrorCode DiagonalDominance 00280 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00281 { 00282 Mat A = (Mat)prob; 00283 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00284 PetscFunctionBegin; 00285 00286 ierr = GetDataID("simple","diagonal-dominance",&id,&has); CHKERRQ(ierr); 00287 HASTOEXIST(has); 00288 ierr = PetscObjectComposedDataGetReal 00289 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00290 if (*flg) rv->r = v; 00291 else { 00292 ierr = compute_dd(A,flg); CHKERRQ(ierr); 00293 if (*flg) { 00294 ierr = PetscObjectComposedDataGetReal 00295 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00296 if (*flg) rv->r = v; 00297 } 00298 } 00299 00300 PetscFunctionReturn(0); 00301 } 00302 00303 /* 00304 * Norms of symmetric / unsymmetric part 00305 */ 00306 #undef __FUNCT__ 00307 #define __FUNCT__ "MatSymmPartNormInf_seq" 00308 /*! 00309 This is where the real work of MatSymmPartNormInf() is done 00310 00311 A few temporary arrays of the 00312 size of a vector are allocated and freed. Otherwise this is just a whole 00313 bunch of pointer chasing. 00314 */ 00315 static PetscErrorCode MatSymmPartNormInf_seq 00316 (Mat A,PetscReal *sn_r,PetscReal *an_r,PetscReal *srn_r,PetscReal *arn_r, 00317 int *nun_r,PetscTruth *done) 00318 { 00319 Mat_SeqAIJ *aij; 00320 int *adx,*bdx,*aii,*bii,*aptr,*bptr,ma,na,i,nun; PetscScalar *va,*vb; 00321 PetscReal snorm,anorm,*sronorm,*aronorm; 00322 PetscErrorCode ierr; 00323 PetscFunctionBegin; 00324 00325 aij = (Mat_SeqAIJ *) A->data; 00326 ierr = MatGetLocalSize(A,&ma,&na); CHKERRQ(ierr); 00327 00328 /* two arrays of pointers into the matrix */ 00329 aii = bii = aij->i; 00330 adx = bdx = aij->j; 00331 va = vb = aij->a; 00332 ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr); 00333 ierr = PetscMalloc(ma*sizeof(int),&bptr); CHKERRQ(ierr); 00334 for (i=0; i<ma; i++) aptr[i] = bptr[i] = aii[i]; 00335 00336 /* temp storage for the inf norm */ 00337 ierr = PetscMalloc(ma*sizeof(PetscReal),&sronorm); CHKERRQ(ierr); 00338 ierr = PetscMalloc(ma*sizeof(PetscReal),&aronorm); CHKERRQ(ierr); 00339 ierr = PetscMemzero(sronorm,ma*sizeof(PetscReal)); CHKERRQ(ierr); 00340 ierr = PetscMemzero(aronorm,ma*sizeof(PetscReal)); CHKERRQ(ierr); 00341 00342 /* update row & frobenius norms with a(i,j)=v and a(j,i)=w */ 00343 #if defined(ANAMODDEBUG) 00344 #define RC_NORM(i,j,v,w) \ 00345 printf("s/a norm calculation: (i,j)=(%d,%d), (vij,vji)=(%e,%e)\n",i,j,v,w); \ 00346 { \ 00347 PetscReal sum = PetscAbsScalar(v+w)/2, dif = PetscAbsScalar(v-w)/2; \ 00348 ASSERT3(!(i==j && v!=w),"same elements on diagonal",i,j); \ 00349 sronorm[i] += sum; if (i!=j) sronorm[j] += sum; \ 00350 aronorm[i] += dif; if (i!=j) aronorm[j] += dif; \ 00351 if (i==j) snorm += sum*sum; \ 00352 else {snorm += 2*sum*sum; anorm += 2*dif*dif;} \ 00353 } 00354 #else 00355 #define RC_NORM(i,j,v,w) \ 00356 { \ 00357 PetscReal sum = PetscAbsScalar(v+w)/2, dif = PetscAbsScalar(v-w)/2; \ 00358 ASSERT3(!(i==j && v!=w),"same elements on diagonal",i,j); \ 00359 sronorm[i] += sum; if (i!=j) sronorm[j] += sum; \ 00360 aronorm[i] += dif; if (i!=j) aronorm[j] += dif; \ 00361 if (i==j) snorm += sum*sum; \ 00362 else {snorm += 2*sum*sum; anorm += 2*dif*dif;} \ 00363 } 00364 #endif 00365 #define ASSERT3(v,s,i,j) \ 00366 if (!(v)) SETERRQ3(1,"Violation of <%s> @ (%d,%d): %d",s,i,j); 00367 #define ASSERT(v,s,i,j,t) \ 00368 if (!(v)) SETERRQ4(1,"Violation of <%s> @ (%d,%d): %d",s,i,j,t); 00369 00370 snorm = anorm = 0; nun = 0; 00371 for (i=0; i<ma; i++) { 00372 int idc,idr; PetscScalar vc,vr; 00373 /* get the a-pointer into the upper triangle */ 00374 while (aptr[i]<aii[i+1] && adx[aptr[i]]<i) aptr[i]++; 00375 ASSERT3(aptr[i]==aii[i+1]||adx[aptr[i]]>=i,"A in upper0",i,adx[aptr[i]]); 00376 /* incorporate L elements that are not structurally symmetric */ 00377 while (bptr[i]<bii[i+1] && bdx[bptr[i]]<i) { 00378 vr = PetscAbsScalar(vb[bptr[i]]); nun++; 00379 RC_NORM(i,bdx[bptr[i]],vr,0.); 00380 bptr[i]++; 00381 } 00382 ASSERT3(bptr[i]==bii[i+1] || bdx[bptr[i]]>=i,"B part done",i,bdx[bptr[i]]); 00383 while (aptr[i]<aii[i+1]) { 00384 /* in row i we are at element j = idc @ adx( aptr[i] ); 00385 vc = A( i, j ) */ 00386 idc = adx[aptr[i]]; vc = va[aptr[i]]; 00387 ASSERT(idc>=i,"A in upper",i,idc,i-idc); 00388 /* in row j=idc look at elements i' = idr @ bdx(bptr[idc]); vr = B(j,i') 00389 first catch up with structurally unsymmetric elements */ 00390 idr = bdx[bptr[idc]]; vr = vb[bptr[idc]]; 00391 while (idr<i) { 00392 RC_NORM(idc,idr,vr,0.); 00393 bptr[idc]++; idr = bdx[bptr[idc]]; vr = vb[bptr[idc]]; nun++; 00394 } 00395 /* now see if (i,j) == (j,i') */ 00396 if (idr==i) { 00397 /* structurally symmetric; process these elements of L and U */ 00398 RC_NORM(i,idc,vc,vr); 00399 bptr[idc]++; 00400 } else { 00401 /* element (j,i') in L corresponds to a later row in U, 00402 so we have a structurally unsymmetric element in U */ 00403 nun++; RC_NORM(i,idc,vc,0.); 00404 } 00405 aptr[i]++; 00406 } 00407 ASSERT(aptr[i]==aii[i+1],"Aptr has traversed row i",i+1,adx[aptr[i]],aii[i+1]); 00408 } 00409 00410 ierr = PetscFree(aptr); CHKERRQ(ierr); 00411 ierr = PetscFree(bptr); CHKERRQ(ierr); 00412 00413 snorm = sqrt(snorm); anorm = sqrt(anorm); 00414 00415 for (i=1; i<ma; i++) { 00416 if (sronorm[i]>sronorm[0]) sronorm[0] = sronorm[i]; 00417 if (aronorm[i]>aronorm[0]) aronorm[0] = aronorm[i]; 00418 } 00419 00420 *sn_r = snorm; *an_r = anorm; *srn_r = sronorm[0]; *arn_r = aronorm[0]; 00421 *nun_r = nun; 00422 *done = PETSC_TRUE; 00423 00424 ierr = PetscFree(sronorm); CHKERRQ(ierr); 00425 ierr = PetscFree(aronorm); CHKERRQ(ierr); 00426 00427 PetscFunctionReturn(0); 00428 } 00429 00430 #undef __FUNCT__ 00431 #define __FUNCT__ "MatSymmPartNormInf" 00432 /*! 00433 This is the computational routine for all quantities relating to 00434 symmetric and anti-symmetric parts of the matrix. 00435 00436 This computation can only be done in single-processor mode. 00437 If the commandline 00438 option (see \ref options) \c "-anamod-force sequential" is specified, 00439 the matrix is gathered on processor zero to compute this quantity. 00440 00441 See MatSymmPartNormInf_seq() for the actual computation. 00442 */ 00443 static PetscErrorCode MatSymmPartNormInf(Mat A,PetscTruth *done) 00444 { 00445 MPI_Comm comm; Mat Awork; int id,nun; 00446 PetscReal snorm,anorm,sronorm,aronorm; 00447 PetscTruth has,alloc_work,do_work,any_work; PetscErrorCode ierr; 00448 00449 PetscFunctionBegin; 00450 00451 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00452 ierr = AnaModGetSequentialMatrix 00453 (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr); 00454 if (!any_work) { 00455 *done = PETSC_FALSE; PetscFunctionReturn(0); 00456 } 00457 if (do_work) { 00458 ierr = MatSymmPartNormInf_seq 00459 (Awork,&snorm,&anorm,&sronorm,&aronorm,&nun,done); CHKERRQ(ierr); 00460 } 00461 if (alloc_work) { 00462 ierr = MatDestroy(Awork); CHKERRQ(ierr); 00463 } 00464 00465 MPI_Bcast((void*)done,1,MPI_INT,0,comm); 00466 MPI_Bcast((void*)&snorm,1,MPIU_REAL,0,comm); 00467 MPI_Bcast((void*)&anorm,1,MPIU_REAL,0,comm); 00468 MPI_Bcast((void*)&sronorm,1,MPIU_REAL,0,comm); 00469 MPI_Bcast((void*)&aronorm,1,MPIU_REAL,0,comm); 00470 MPI_Bcast((void*)&nun,1,MPI_INT,0,comm); 00471 00472 /* infinity norms of symmetric / anti-symmetric part */ 00473 ierr = GetDataID("simple","symmetry-snorm",&id,&has); CHKERRQ(ierr); 00474 HASTOEXIST(has); 00475 ierr = PetscObjectComposedDataSetReal 00476 ((PetscObject)A,id,sronorm); CHKERRQ(ierr); 00477 00478 ierr = GetDataID("simple","symmetry-anorm",&id,&has); CHKERRQ(ierr); 00479 HASTOEXIST(has); 00480 ierr = PetscObjectComposedDataSetReal 00481 ((PetscObject)A,id,aronorm); CHKERRQ(ierr); 00482 00483 /* frobenius norms of symmetric / anti-symmetric part */ 00484 ierr = GetDataID("simple","symmetry-fsnorm",&id,&has); CHKERRQ(ierr); 00485 HASTOEXIST(has); 00486 ierr = PetscObjectComposedDataSetReal 00487 ((PetscObject)A,id,snorm); CHKERRQ(ierr); 00488 00489 ierr = GetDataID("simple","symmetry-fanorm",&id,&has); CHKERRQ(ierr); 00490 HASTOEXIST(has); 00491 ierr = PetscObjectComposedDataSetReal 00492 ((PetscObject)A,id,anorm); CHKERRQ(ierr); 00493 00494 ierr = GetDataID("structure","n-struct-unsymm",&id,&has); CHKERRQ(ierr); 00495 HASTOEXIST(has); 00496 ierr = PetscObjectComposedDataSetInt 00497 ((PetscObject)A,id,nun); CHKERRQ(ierr); 00498 00499 PetscFunctionReturn(0); 00500 } 00501 00502 #undef __FUNCT__ 00503 #define __FUNCT__ "SymmetrySNorm" 00504 /*! Compute and store the infinity norm of the symmetric part of the matrix. 00505 The actual computation is done in MatSymmPartNormInf(); 00506 see that routine for remarks on parallelism. 00507 */ 00508 static PetscErrorCode SymmetrySNorm 00509 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00510 { 00511 Mat A = (Mat)prob; 00512 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00513 PetscFunctionBegin; 00514 ierr = GetDataID("simple","symmetry-snorm",&id,&has); CHKERRQ(ierr); 00515 HASTOEXIST(has); 00516 ierr = PetscObjectComposedDataGetReal 00517 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00518 if (*flg) rv->r = v; 00519 else { 00520 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00521 if (*flg) { 00522 ierr = PetscObjectComposedDataGetReal 00523 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00524 if (*flg) rv->r = v; 00525 } 00526 } 00527 00528 PetscFunctionReturn(0); 00529 } 00530 00531 #undef __FUNCT__ 00532 #define __FUNCT__ "SymmetryANorm" 00533 /*! Compute and store the infinity norm of the antisymmetric part of the matrix. 00534 The actual computation is done in MatSymmPartNormInf(); 00535 see that routine for remarks on parallelism. 00536 */ 00537 static PetscErrorCode SymmetryANorm 00538 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00539 { 00540 Mat A = (Mat)prob; 00541 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00542 PetscFunctionBegin; 00543 ierr = GetDataID("simple","symmetry-anorm",&id,&has); CHKERRQ(ierr); 00544 HASTOEXIST(has); 00545 ierr = PetscObjectComposedDataGetReal 00546 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00547 if (*flg) rv->r = v; 00548 else { 00549 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00550 if (*flg) { 00551 ierr = PetscObjectComposedDataGetReal 00552 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00553 if (*flg) rv->r = v; 00554 } 00555 } 00556 00557 PetscFunctionReturn(0); 00558 } 00559 00560 #undef __FUNCT__ 00561 #define __FUNCT__ "SymmetryFSNorm" 00562 /*! Compute and store the Frobenius norm of the symmetric part of the matrix. 00563 The actual computation is done in MatSymmPartNormInf(); 00564 see that routine for remarks on parallelism. 00565 */ 00566 static PetscErrorCode SymmetryFSNorm 00567 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00568 { 00569 Mat A = (Mat)prob; 00570 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00571 PetscFunctionBegin; 00572 ierr = GetDataID("simple","symmetry-fsnorm",&id,&has); CHKERRQ(ierr); 00573 HASTOEXIST(has); 00574 ierr = PetscObjectComposedDataGetReal 00575 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00576 if (*flg) rv->r = v; 00577 else { 00578 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00579 if (*flg) { 00580 ierr = PetscObjectComposedDataGetReal 00581 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00582 if (*flg) rv->r = v; 00583 } 00584 } 00585 00586 PetscFunctionReturn(0); 00587 } 00588 00589 #undef __FUNCT__ 00590 #define __FUNCT__ "SymmetryFANorm" 00591 /*! Compute and store the Frobenius norm of the antisymmetric part of the matrix. 00592 The actual computation is done in MatSymmPartNormInf(); 00593 see that routine for remarks on parallelism. 00594 */ 00595 static PetscErrorCode SymmetryFANorm 00596 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00597 { 00598 Mat A = (Mat)prob; 00599 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00600 PetscFunctionBegin; 00601 ierr = GetDataID("simple","symmetry-fanorm",&id,&has); CHKERRQ(ierr); 00602 HASTOEXIST(has); 00603 ierr = PetscObjectComposedDataGetReal 00604 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00605 if (*flg) rv->r = v; 00606 else { 00607 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00608 if (*flg) { 00609 ierr = PetscObjectComposedDataGetReal 00610 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00611 if (*flg) rv->r = v; 00612 } 00613 } 00614 00615 PetscFunctionReturn(0); 00616 } 00617 00618 #undef __FUNCT__ 00619 #define __FUNCT__ "NUnstruct" 00620 /*! Compute and store the number of structurally unsymmetric elements of a matrix. 00621 The actual computation is done in MatSymmPartNormInf(); 00622 see that routine for remarks on parallelism. 00623 00624 The "n-struct-unsymm" feature is declared to be in category "structure", 00625 rather than simple. 00626 */ 00627 static PetscErrorCode NUnstruct 00628 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00629 { 00630 Mat A = (Mat)prob; 00631 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00632 PetscFunctionBegin; 00633 ierr = GetDataID("structure","n-struct-unsymm",&id,&has); CHKERRQ(ierr); 00634 HASTOEXIST(has); 00635 ierr = PetscObjectComposedDataGetInt 00636 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00637 if (*flg) rv->i = v; 00638 else { 00639 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00640 if (*flg) { 00641 ierr = PetscObjectComposedDataGetInt 00642 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00643 if (*flg) rv->i = v; 00644 } 00645 } 00646 00647 PetscFunctionReturn(0); 00648 } 00649 00650 #undef __FUNCT__ 00651 #define __FUNCT__ "RegisterSimpleModules" 00652 /*! 00653 Declare norm-like modules that can be performed with simple calculations. 00654 */ 00655 PetscErrorCode RegisterSimpleModules() 00656 { 00657 PetscErrorCode ierr; 00658 PetscFunctionBegin; 00659 00660 ierr = RegisterModule 00661 ("simple","trace",ANALYSISDOUBLE,&Trace); CHKERRQ(ierr); 00662 ierr = RegisterModule 00663 ("simple","trace-abs",ANALYSISDOUBLE,&TraceAbs); CHKERRQ(ierr); 00664 00665 ierr = RegisterModule 00666 ("simple","norm1",ANALYSISDOUBLE,&norm1); CHKERRQ(ierr); 00667 ierr = RegisterModule 00668 ("simple","normInf",ANALYSISDOUBLE,&normInf); CHKERRQ(ierr); 00669 ierr = RegisterModule 00670 ("simple","normF",ANALYSISDOUBLE,&normF); CHKERRQ(ierr); 00671 00672 ierr = RegisterModule 00673 ("simple","diagonal-dominance",ANALYSISDOUBLE,&DiagonalDominance); CHKERRQ(ierr); 00674 00675 ierr = RegisterModule 00676 ("simple","symmetry-snorm",ANALYSISDOUBLE,&SymmetrySNorm); CHKERRQ(ierr); 00677 ierr = RegisterModule 00678 ("simple","symmetry-anorm",ANALYSISDOUBLE,&SymmetryANorm); CHKERRQ(ierr); 00679 ierr = RegisterModule 00680 ("simple","symmetry-fsnorm",ANALYSISDOUBLE,&SymmetryFSNorm); CHKERRQ(ierr); 00681 ierr = RegisterModule 00682 ("simple","symmetry-fanorm",ANALYSISDOUBLE,&SymmetryFANorm); CHKERRQ(ierr); 00683 00684 ierr = RegisterModule 00685 ("structure","n-struct-unsymm",ANALYSISINTEGER,&NUnstruct); CHKERRQ(ierr); 00686 00687 PetscFunctionReturn(0); 00688 } 00689 00690 #undef __FUNCT__ 00691 #define __FUNCT__ "DeRegisterSimpleModules" 00692 PetscErrorCode DeRegisterSimpleModules(void) 00693 { 00694 PetscFunctionBegin; 00695 PetscFunctionReturn(0); 00696 }