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 PetscBool 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,PetscBool *flg) 00107 { 00108 Mat A = (Mat)prob; 00109 PetscReal v=0; PetscBool 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(MPI_COMM_WORLD,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,PetscBool *flg) 00136 { 00137 Mat A = (Mat)prob; 00138 PetscReal v=0; PetscBool 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,PetscBool *flg) 00160 { 00161 Mat A = (Mat)prob; 00162 PetscReal v=0; PetscBool 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,PetscBool *flg) 00187 { 00188 Mat A = (Mat)prob; 00189 PetscReal v=0; PetscBool 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,PetscBool *flg) 00214 { 00215 Mat A = (Mat)prob; 00216 PetscReal v=0; PetscBool 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,PetscBool *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 ierr = MPI_Allreduce((void*)&r,(void*)&gr,1,MPIU_REAL,MPI_MIN,comm); 00268 { 00269 int id; PetscBool 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 #undef __FUNCT__ 00280 #define __FUNCT__ "DiagonalDominance" 00281 static PetscErrorCode DiagonalDominance 00282 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00283 { 00284 Mat A = (Mat)prob; 00285 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00286 PetscFunctionBegin; 00287 00288 ierr = GetDataID("simple","diagonal-dominance",&id,&has); CHKERRQ(ierr); 00289 HASTOEXIST(has); 00290 ierr = PetscObjectComposedDataGetReal 00291 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00292 if (*flg) rv->r = v; 00293 else { 00294 ierr = compute_dd(A,flg); CHKERRQ(ierr); 00295 if (*flg) { 00296 ierr = PetscObjectComposedDataGetReal 00297 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00298 if (*flg) rv->r = v; 00299 } 00300 } 00301 00302 PetscFunctionReturn(0); 00303 } 00304 00305 /* 00306 * Norms of symmetric / unsymmetric part 00307 */ 00308 #undef __FUNCT__ 00309 #define __FUNCT__ "MatSymmPartNormInf_seq" 00310 /*! 00311 This is where the real work of MatSymmPartNormInf() is done 00312 00313 A few temporary arrays of the 00314 size of a vector are allocated and freed. Otherwise this is just a whole 00315 bunch of pointer chasing. 00316 */ 00317 static PetscErrorCode MatSymmPartNormInf_seq 00318 (Mat A,PetscReal *sn_r,PetscReal *an_r,PetscReal *srn_r,PetscReal *arn_r, 00319 int *nun_r,PetscBool *done) 00320 { 00321 Mat_SeqAIJ *aij; 00322 int *adx,*bdx,*aii,*bii,*aptr,*bptr,ma,na,i,nun; PetscScalar *va,*vb; 00323 PetscReal snorm,anorm,*sronorm,*aronorm; 00324 PetscErrorCode ierr; 00325 PetscFunctionBegin; 00326 00327 aij = (Mat_SeqAIJ *) A->data; 00328 ierr = MatGetLocalSize(A,&ma,&na); CHKERRQ(ierr); 00329 00330 /* two arrays of pointers into the matrix */ 00331 aii = bii = aij->i; 00332 adx = bdx = aij->j; 00333 va = vb = aij->a; 00334 ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr); 00335 ierr = PetscMalloc(ma*sizeof(int),&bptr); CHKERRQ(ierr); 00336 for (i=0; i<ma; i++) aptr[i] = bptr[i] = aii[i]; 00337 00338 /* temp storage for the inf norm */ 00339 ierr = PetscMalloc(ma*sizeof(PetscReal),&sronorm); CHKERRQ(ierr); 00340 ierr = PetscMalloc(ma*sizeof(PetscReal),&aronorm); CHKERRQ(ierr); 00341 ierr = PetscMemzero(sronorm,ma*sizeof(PetscReal)); CHKERRQ(ierr); 00342 ierr = PetscMemzero(aronorm,ma*sizeof(PetscReal)); CHKERRQ(ierr); 00343 00344 /* update row & frobenius norms with a(i,j)=v and a(j,i)=w */ 00345 #if defined(ANAMODDEBUG) 00346 #define RC_NORM(i,j,v,w) \ 00347 printf("s/a norm calculation: (i,j)=(%d,%d), (vij,vji)=(%e,%e)\n",i,j,v,w); \ 00348 { \ 00349 PetscReal sum = PetscAbsScalar(v+w)/2, dif = PetscAbsScalar(v-w)/2; \ 00350 ASSERT3(!(i==j && v!=w),"same elements on diagonal",i,j); \ 00351 sronorm[i] += sum; if (i!=j) sronorm[j] += sum; \ 00352 aronorm[i] += dif; if (i!=j) aronorm[j] += dif; \ 00353 if (i==j) snorm += sum*sum; \ 00354 else {snorm += 2*sum*sum; anorm += 2*dif*dif;} \ 00355 } 00356 #else 00357 #define RC_NORM(i,j,v,w) \ 00358 { \ 00359 PetscReal sum = PetscAbsScalar(v+w)/2, dif = PetscAbsScalar(v-w)/2; \ 00360 ASSERT3(!(i==j && v!=w),"same elements on diagonal",i,j); \ 00361 sronorm[i] += sum; if (i!=j) sronorm[j] += sum; \ 00362 aronorm[i] += dif; if (i!=j) aronorm[j] += dif; \ 00363 if (i==j) snorm += sum*sum; \ 00364 else {snorm += 2*sum*sum; anorm += 2*dif*dif;} \ 00365 } 00366 #endif 00367 #define ASSERT3(v,s,i,j) \ 00368 if (!(v)) SETERRQ3(MPI_COMM_WORLD,1,"Violation of <%s> @ (%d,%d): %d",s,i,j); 00369 #define ASSERT(v,s,i,j,t) \ 00370 if (!(v)) SETERRQ4(MPI_COMM_WORLD,1,"Violation of <%s> @ (%d,%d): %d",s,i,j,t); 00371 00372 snorm = anorm = 0; nun = 0; 00373 for (i=0; i<ma; i++) { 00374 int idc,idr; PetscScalar vc,vr; 00375 /* get the a-pointer into the upper triangle */ 00376 while (aptr[i]<aii[i+1] && adx[aptr[i]]<i) aptr[i]++; 00377 ASSERT3(aptr[i]==aii[i+1]||adx[aptr[i]]>=i,"A in upper0",i,adx[aptr[i]]); 00378 /* incorporate L elements that are not structurally symmetric */ 00379 while (bptr[i]<bii[i+1] && bdx[bptr[i]]<i) { 00380 vr = PetscAbsScalar(vb[bptr[i]]); nun++; 00381 RC_NORM(i,bdx[bptr[i]],vr,0.); 00382 bptr[i]++; 00383 } 00384 ASSERT3(bptr[i]==bii[i+1] || bdx[bptr[i]]>=i,"B part done",i,bdx[bptr[i]]); 00385 while (aptr[i]<aii[i+1]) { 00386 /* in row i we are at element j = idc @ adx( aptr[i] ); 00387 vc = A( i, j ) */ 00388 idc = adx[aptr[i]]; vc = va[aptr[i]]; 00389 ASSERT(idc>=i,"A in upper",i,idc,i-idc); 00390 /* in row j=idc look at elements i' = idr @ bdx(bptr[idc]); vr = B(j,i') 00391 first catch up with structurally unsymmetric elements */ 00392 idr = bdx[bptr[idc]]; vr = vb[bptr[idc]]; 00393 while (idr<i) { 00394 RC_NORM(idc,idr,vr,0.); 00395 bptr[idc]++; idr = bdx[bptr[idc]]; vr = vb[bptr[idc]]; nun++; 00396 } 00397 /* now see if (i,j) == (j,i') */ 00398 if (idr==i) { 00399 /* structurally symmetric; process these elements of L and U */ 00400 RC_NORM(i,idc,vc,vr); 00401 bptr[idc]++; 00402 } else { 00403 /* element (j,i') in L corresponds to a later row in U, 00404 so we have a structurally unsymmetric element in U */ 00405 nun++; RC_NORM(i,idc,vc,0.); 00406 } 00407 aptr[i]++; 00408 } 00409 ASSERT(aptr[i]==aii[i+1],"Aptr has traversed row i",i+1,adx[aptr[i]],aii[i+1]); 00410 } 00411 00412 ierr = PetscFree(aptr); CHKERRQ(ierr); 00413 ierr = PetscFree(bptr); CHKERRQ(ierr); 00414 00415 snorm = sqrt(snorm); anorm = sqrt(anorm); 00416 00417 for (i=1; i<ma; i++) { 00418 if (sronorm[i]>sronorm[0]) sronorm[0] = sronorm[i]; 00419 if (aronorm[i]>aronorm[0]) aronorm[0] = aronorm[i]; 00420 } 00421 00422 *sn_r = snorm; *an_r = anorm; *srn_r = sronorm[0]; *arn_r = aronorm[0]; 00423 *nun_r = nun; 00424 *done = PETSC_TRUE; 00425 00426 ierr = PetscFree(sronorm); CHKERRQ(ierr); 00427 ierr = PetscFree(aronorm); CHKERRQ(ierr); 00428 00429 PetscFunctionReturn(0); 00430 } 00431 00432 #undef __FUNCT__ 00433 #define __FUNCT__ "MatSymmPartNormInf" 00434 /*! 00435 This is the computational routine for all quantities relating to 00436 symmetric and anti-symmetric parts of the matrix. 00437 00438 This computation can only be done in single-processor mode. 00439 If the commandline 00440 option (see \ref options) \c "-anamod-force sequential" is specified, 00441 the matrix is gathered on processor zero to compute this quantity. 00442 00443 See MatSymmPartNormInf_seq() for the actual computation. 00444 */ 00445 static PetscErrorCode MatSymmPartNormInf(Mat A,PetscBool *done) 00446 { 00447 MPI_Comm comm; Mat Awork; int id,nun; 00448 PetscReal snorm,anorm,sronorm,aronorm; 00449 PetscBool has,alloc_work,do_work,any_work; PetscErrorCode ierr; 00450 00451 PetscFunctionBegin; 00452 00453 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00454 ierr = AnaModGetSequentialMatrix 00455 (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr); 00456 if (!any_work) { 00457 *done = PETSC_FALSE; PetscFunctionReturn(0); 00458 } 00459 if (do_work) { 00460 ierr = MatSymmPartNormInf_seq 00461 (Awork,&snorm,&anorm,&sronorm,&aronorm,&nun,done); CHKERRQ(ierr); 00462 } 00463 if (alloc_work) { 00464 ierr = MatDestroy(&Awork); CHKERRQ(ierr); 00465 } 00466 00467 ierr = MPI_Bcast((void*)done,1,MPI_INT,0,comm); 00468 ierr = MPI_Bcast((void*)&snorm,1,MPIU_REAL,0,comm); 00469 ierr = MPI_Bcast((void*)&anorm,1,MPIU_REAL,0,comm); 00470 ierr = MPI_Bcast((void*)&sronorm,1,MPIU_REAL,0,comm); 00471 ierr = MPI_Bcast((void*)&aronorm,1,MPIU_REAL,0,comm); 00472 ierr = MPI_Bcast((void*)&nun,1,MPI_INT,0,comm); 00473 00474 /* infinity norms of symmetric / anti-symmetric part */ 00475 ierr = GetDataID("simple","symmetry-snorm",&id,&has); CHKERRQ(ierr); 00476 HASTOEXIST(has); 00477 ierr = PetscObjectComposedDataSetReal 00478 ((PetscObject)A,id,sronorm); CHKERRQ(ierr); 00479 00480 ierr = GetDataID("simple","symmetry-anorm",&id,&has); CHKERRQ(ierr); 00481 HASTOEXIST(has); 00482 ierr = PetscObjectComposedDataSetReal 00483 ((PetscObject)A,id,aronorm); CHKERRQ(ierr); 00484 00485 /* frobenius norms of symmetric / anti-symmetric part */ 00486 ierr = GetDataID("simple","symmetry-fsnorm",&id,&has); CHKERRQ(ierr); 00487 HASTOEXIST(has); 00488 ierr = PetscObjectComposedDataSetReal 00489 ((PetscObject)A,id,snorm); CHKERRQ(ierr); 00490 00491 ierr = GetDataID("simple","symmetry-fanorm",&id,&has); CHKERRQ(ierr); 00492 HASTOEXIST(has); 00493 ierr = PetscObjectComposedDataSetReal 00494 ((PetscObject)A,id,anorm); CHKERRQ(ierr); 00495 00496 ierr = GetDataID("structure","n-struct-unsymm",&id,&has); CHKERRQ(ierr); 00497 HASTOEXIST(has); 00498 ierr = PetscObjectComposedDataSetInt 00499 ((PetscObject)A,id,nun); CHKERRQ(ierr); 00500 00501 PetscFunctionReturn(0); 00502 } 00503 00504 #undef __FUNCT__ 00505 #define __FUNCT__ "SymmetrySNorm" 00506 /*! Compute and store the infinity norm of the symmetric part of the matrix. 00507 The actual computation is done in MatSymmPartNormInf(); 00508 see that routine for remarks on parallelism. 00509 */ 00510 static PetscErrorCode SymmetrySNorm 00511 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00512 { 00513 Mat A = (Mat)prob; 00514 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00515 PetscFunctionBegin; 00516 ierr = GetDataID("simple","symmetry-snorm",&id,&has); CHKERRQ(ierr); 00517 HASTOEXIST(has); 00518 ierr = PetscObjectComposedDataGetReal 00519 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00520 if (*flg) rv->r = v; 00521 else { 00522 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00523 if (*flg) { 00524 ierr = PetscObjectComposedDataGetReal 00525 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00526 if (*flg) rv->r = v; 00527 } 00528 } 00529 00530 PetscFunctionReturn(0); 00531 } 00532 00533 #undef __FUNCT__ 00534 #define __FUNCT__ "SymmetryANorm" 00535 /*! Compute and store the infinity norm of the antisymmetric part of the matrix. 00536 The actual computation is done in MatSymmPartNormInf(); 00537 see that routine for remarks on parallelism. 00538 */ 00539 static PetscErrorCode SymmetryANorm 00540 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00541 { 00542 Mat A = (Mat)prob; 00543 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00544 PetscFunctionBegin; 00545 ierr = GetDataID("simple","symmetry-anorm",&id,&has); CHKERRQ(ierr); 00546 HASTOEXIST(has); 00547 ierr = PetscObjectComposedDataGetReal 00548 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00549 if (*flg) rv->r = v; 00550 else { 00551 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00552 if (*flg) { 00553 ierr = PetscObjectComposedDataGetReal 00554 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00555 if (*flg) rv->r = v; 00556 } 00557 } 00558 00559 PetscFunctionReturn(0); 00560 } 00561 00562 #undef __FUNCT__ 00563 #define __FUNCT__ "SymmetryFSNorm" 00564 /*! Compute and store the Frobenius norm of the symmetric part of the matrix. 00565 The actual computation is done in MatSymmPartNormInf(); 00566 see that routine for remarks on parallelism. 00567 */ 00568 static PetscErrorCode SymmetryFSNorm 00569 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00570 { 00571 Mat A = (Mat)prob; 00572 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00573 PetscFunctionBegin; 00574 ierr = GetDataID("simple","symmetry-fsnorm",&id,&has); CHKERRQ(ierr); 00575 HASTOEXIST(has); 00576 ierr = PetscObjectComposedDataGetReal 00577 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00578 if (*flg) rv->r = v; 00579 else { 00580 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00581 if (*flg) { 00582 ierr = PetscObjectComposedDataGetReal 00583 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00584 if (*flg) rv->r = v; 00585 } 00586 } 00587 00588 PetscFunctionReturn(0); 00589 } 00590 00591 #undef __FUNCT__ 00592 #define __FUNCT__ "SymmetryFANorm" 00593 /*! Compute and store the Frobenius norm of the antisymmetric part of the matrix. 00594 The actual computation is done in MatSymmPartNormInf(); 00595 see that routine for remarks on parallelism. 00596 */ 00597 static PetscErrorCode SymmetryFANorm 00598 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00599 { 00600 Mat A = (Mat)prob; 00601 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00602 PetscFunctionBegin; 00603 ierr = GetDataID("simple","symmetry-fanorm",&id,&has); CHKERRQ(ierr); 00604 HASTOEXIST(has); 00605 ierr = PetscObjectComposedDataGetReal 00606 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00607 if (*flg) rv->r = v; 00608 else { 00609 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00610 if (*flg) { 00611 ierr = PetscObjectComposedDataGetReal 00612 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00613 if (*flg) rv->r = v; 00614 } 00615 } 00616 00617 PetscFunctionReturn(0); 00618 } 00619 00620 #undef __FUNCT__ 00621 #define __FUNCT__ "NUnstruct" 00622 /*! Compute and store the number of structurally unsymmetric elements of a matrix. 00623 The actual computation is done in MatSymmPartNormInf(); 00624 see that routine for remarks on parallelism. 00625 00626 The "n-struct-unsymm" feature is declared to be in category "structure", 00627 rather than simple. 00628 */ 00629 static PetscErrorCode NUnstruct 00630 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00631 { 00632 Mat A = (Mat)prob; 00633 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00634 PetscFunctionBegin; 00635 ierr = GetDataID("structure","n-struct-unsymm",&id,&has); CHKERRQ(ierr); 00636 HASTOEXIST(has); 00637 ierr = PetscObjectComposedDataGetInt 00638 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00639 if (*flg) rv->i = v; 00640 else { 00641 ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr); 00642 if (*flg) { 00643 ierr = PetscObjectComposedDataGetInt 00644 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00645 if (*flg) rv->i = v; 00646 } 00647 } 00648 00649 PetscFunctionReturn(0); 00650 } 00651 00652 #undef __FUNCT__ 00653 #define __FUNCT__ "RegisterSimpleModules" 00654 /*! 00655 Declare norm-like modules that can be performed with simple calculations. 00656 */ 00657 PetscErrorCode RegisterSimpleModules() 00658 { 00659 PetscErrorCode ierr; 00660 PetscFunctionBegin; 00661 00662 ierr = RegisterModule 00663 ("simple","trace",ANALYSISDOUBLE,&Trace); CHKERRQ(ierr); 00664 ierr = RegisterModule 00665 ("simple","trace-abs",ANALYSISDOUBLE,&TraceAbs); CHKERRQ(ierr); 00666 00667 ierr = RegisterModule 00668 ("simple","norm1",ANALYSISDOUBLE,&norm1); CHKERRQ(ierr); 00669 ierr = RegisterModule 00670 ("simple","normInf",ANALYSISDOUBLE,&normInf); CHKERRQ(ierr); 00671 ierr = RegisterModule 00672 ("simple","normF",ANALYSISDOUBLE,&normF); CHKERRQ(ierr); 00673 00674 ierr = RegisterModule 00675 ("simple","diagonal-dominance",ANALYSISDOUBLE,&DiagonalDominance); CHKERRQ(ierr); 00676 00677 ierr = RegisterModule 00678 ("simple","symmetry-snorm",ANALYSISDOUBLE,&SymmetrySNorm); CHKERRQ(ierr); 00679 ierr = RegisterModule 00680 ("simple","symmetry-anorm",ANALYSISDOUBLE,&SymmetryANorm); CHKERRQ(ierr); 00681 ierr = RegisterModule 00682 ("simple","symmetry-fsnorm",ANALYSISDOUBLE,&SymmetryFSNorm); CHKERRQ(ierr); 00683 ierr = RegisterModule 00684 ("simple","symmetry-fanorm",ANALYSISDOUBLE,&SymmetryFANorm); CHKERRQ(ierr); 00685 00686 ierr = RegisterModule 00687 ("structure","n-struct-unsymm",ANALYSISINTEGER,&NUnstruct); CHKERRQ(ierr); 00688 00689 PetscFunctionReturn(0); 00690 } 00691 00692 #undef __FUNCT__ 00693 #define __FUNCT__ "DeRegisterSimpleModules" 00694 PetscErrorCode DeRegisterSimpleModules(void) 00695 { 00696 PetscFunctionBegin; 00697 PetscFunctionReturn(0); 00698 }