MITgcm |
|
|
|
|||
File indexing completed on 2023-11-05 05:10:01 UTC
view on githubraw file Latest commit 65754df4 on 2023-11-04 17:55:24 UTCb35bd3101a Jean*0001 #include "ADMTLM_OPTIONS.h" a5276edbc9 Patr*0002 ee5ca1ad15 Patr*0003 subroutine admtlm_dsvd ( mythid ) 0004 c b35bd3101a Jean*0005 c This example program is intended to illustrate the ee5ca1ad15 Patr*0006 c the use of ARPACK to compute the Singular Value Decomposition. b35bd3101a Jean*0007 c ee5ca1ad15 Patr*0008 c This code shows how to use ARPACK to find a few of the b35bd3101a Jean*0009 c largest singular values(sigma) and corresponding right singular ee5ca1ad15 Patr*0010 c vectors (v) for the the matrix A by solving the symmetric problem: b35bd3101a Jean*0011 c ee5ca1ad15 Patr*0012 c (A'*A)*v = sigma*v b35bd3101a Jean*0013 c ee5ca1ad15 Patr*0014 c where A is an m by n real matrix. 0015 c 0016 c This code may be easily modified to estimate the 2-norm 0017 c condition number largest(sigma)/smallest(sigma) by setting 0018 c which = 'BE' below. This will ask for a few of the smallest 0019 c and a few of the largest singular values simultaneously. 0020 c The condition number could then be estimated by taking 0021 c the ratio of the largest and smallest singular values. 0022 c 0023 c This formulation is appropriate when m .ge. n. 0024 c Reverse the roles of A and A' in the case that m .le. n. 0025 c b35bd3101a Jean*0026 c The main points illustrated here are ee5ca1ad15 Patr*0027 c b35bd3101a Jean*0028 c 1) How to declare sufficient memory to find NEV 0029 c largest singular values of A . ee5ca1ad15 Patr*0030 c b35bd3101a Jean*0031 c 2) Illustration of the reverse communication interface 0032 c needed to utilize the top level ARPACK routine DSAUPD ee5ca1ad15 Patr*0033 c that computes the quantities needed to construct 0034 c the desired singular values and vectors(if requested). 0035 c 0036 c 3) How to extract the desired singular values and vectors 0037 c using the ARPACK routine DSEUPD. 0038 c b35bd3101a Jean*0039 c 4) How to construct the left singular vectors U from the ee5ca1ad15 Patr*0040 c right singular vectors V to obtain the decomposition 0041 c 0042 c A*V = U*S 0043 c 0044 c where S = diag(sigma_1, sigma_2, ..., sigma_k). 0045 c 0046 c The only thing that must be supplied in order to use this b35bd3101a Jean*0047 c routine on your problem is to change the array dimensions 0048 c appropriately, to specify WHICH singular values you want to 0049 c compute and to supply a the matrix-vector products ee5ca1ad15 Patr*0050 c 0051 c w <- Ax 0052 c y <- A'w 0053 c b35bd3101a Jean*0054 c in place of the calls to AV( ) and ATV( ) respectively below. ee5ca1ad15 Patr*0055 c 0056 c Further documentation is available in the header of DSAUPD 0057 c which may be found in the SRC directory. 0058 c 0059 c This codes implements 0060 c 0061 c\Example-1 0062 c ... Suppose we want to solve A'A*v = sigma*v in regular mode, b35bd3101a Jean*0063 c where A is derived from the simplest finite difference ee5ca1ad15 Patr*0064 c discretization of the 2-dimensional kernel K(s,t)dt where 0065 c 0066 c K(s,t) = s(t-1) if 0 .le. s .le. t .le. 1, b35bd3101a Jean*0067 c t(s-1) if 0 .le. t .lt. s .le. 1. ee5ca1ad15 Patr*0068 c 0069 c See subroutines AV and ATV for details. 0070 c ... OP = A'*A and B = I. 0071 c ... Assume "call av (n,x,y)" computes y = A*x 0072 c ... Assume "call atv (n,y,w)" computes w = A'*y 0073 c ... Assume exact shifts are used 0074 c ... 0075 c 0076 c\BeginLib 0077 c 0078 c\Routines called: 0079 c dsaupd ARPACK reverse communication interface routine. 0080 c dseupd ARPACK routine that returns Ritz values and (optionally) 0081 c Ritz vectors. 0082 c dnrm2 Level 1 BLAS that computes the norm of a vector. 0083 c daxpy Level 1 BLAS that computes y <- alpha*x+y. 0084 c dscal Level 1 BLAS thst computes x <- x*alpha. 0085 c dcopy Level 1 BLAS thst computes y <- x. 0086 c 0087 c\Author 0088 c Richard Lehoucq 0089 c Danny Sorensen 0090 c Chao Yang 0091 c Dept. of Computational & 0092 c Applied Mathematics 0093 c Rice University 0094 c Houston, Texas 0095 c 0096 c\SCCS Information: @(#) 0097 c FILE: svd.F SID: 2.4 DATE OF SID: 10/17/00 RELEASE: 2 0098 c 0099 c\Remarks 0100 c 1. None 0101 c 0102 c\EndLib 0103 c 0104 c----------------------------------------------------------------------- 0105 c 0106 c %------------------------------------------------------% 0107 c | Storage Declarations: | 0108 c | | 0109 c | It is assumed that A is M by N with M .ge. N. | 0110 c | | 0111 c | The maximum dimensions for all arrays are | 0112 c | set here to accommodate a problem size of | 0113 c | M .le. MAXM and N .le. MAXN | 0114 c | | 0115 c | The NEV right singular vectors will be computed in | 0116 c | the N by NCV array V. | 0117 c | | 0118 c | The NEV left singular vectors will be computed in | 0119 c | the M by NEV array U. | 0120 c | | 0121 c | NEV is the number of singular values requested. | 0122 c | See specifications for ARPACK usage below. | 0123 c | | 0124 c | NCV is the largest number of basis vectors that will | 0125 c | be used in the Implicitly Restarted Arnoldi | 0126 c | Process. Work per major iteration is | 0127 c | proportional to N*NCV*NCV. | 0128 c | | 0129 c | You must set: | 0130 c | | 0131 c | MAXM: Maximum number of rows of the A allowed. | 0132 c | MAXN: Maximum number of columns of the A allowed. | 0133 c | MAXNEV: Maximum NEV allowed | 0134 c | MAXNCV: Maximum NCV allowed | 0135 c %------------------------------------------------------% 0136 c a5276edbc9 Patr*0137 C !USES: 0138 C == Global variables === 0139 #include "SIZE.h" 0140 #include "EEPARAMS.h" 0141 #include "PARAMS.h" 0142 #ifdef ALLOW_ADMTLM 7c50f07931 Mart*0143 # ifdef ALLOW_AUTODIFF_TAMC 0144 # include "tamc.h" 0145 # endif 4d72283393 Mart*0146 # include "CTRL.h" 65754df434 Mart*0147 # include "OPTIMCYCLE.h" e4939baa16 Patr*0148 # include "cost.h" 0149 # include "adcost.h" 0150 # include "g_cost.h" a5276edbc9 Patr*0151 #endif 0152 0153 C !INPUT/OUTPUT PARAMETERS: 0154 C == Routine arguments == 0155 INTEGER mythid 0156 0157 C == PARAMETERS == faf44775ba Patr*0158 integer maxnev, maxncv, ldv, ldu ee5ca1ad15 Patr*0159 cph( faf44775ba Patr*0160 parameter ( maxnev=15, maxncv=30, ldu = maxm, ldv=maxn ) a5276edbc9 Patr*0161 character*(80) fnameGlobal ee5ca1ad15 Patr*0162 cph) 0163 c 0164 c %--------------% 0165 c | Local Arrays | 0166 c %--------------% 0167 c 0168 Double precision b35bd3101a Jean*0169 & v(ldv,maxncv), u(ldu, maxnev), 0170 & workl(maxncv*(maxncv+8)), workd(3*maxn), ee5ca1ad15 Patr*0171 & s(maxncv,2), resid(maxn), ax(maxm) 0172 logical select(maxncv) 0173 integer iparam(11), ipntr(11) 0174 c 0175 c %---------------% 0176 c | Local Scalars | 0177 c %---------------% 0178 c 0179 character bmat*1, which*2 0180 integer ido, m, n, nev, ncv, lworkl, info, ierr, 0181 & j, ishfts, maxitr, mode, nconv 0182 logical rvec b35bd3101a Jean*0183 Double precision ee5ca1ad15 Patr*0184 & tol, sigma, temp 0185 c 0186 c %------------% 0187 c | Parameters | 0188 c %------------% 0189 c 0190 Double precision 0191 & one, zero 0192 parameter (one = 1.0D+0, zero = 0.0D+0) b35bd3101a Jean*0193 c ee5ca1ad15 Patr*0194 c %-----------------------------% 0195 c | BLAS & LAPACK routines used | 0196 c %-----------------------------% 0197 c b35bd3101a Jean*0198 Double precision ee5ca1ad15 Patr*0199 & dnrm2 0200 external dnrm2, daxpy, dcopy, dscal 0201 0202 cph( a5276edbc9 Patr*0203 integer l, ilinsysinfo, ii, jj ee5ca1ad15 Patr*0204 integer ipiv(maxn) 0205 integer phiwork(maxn) 0206 0207 double precision ferr, berr 0208 double precision phtmpin(maxn), phtmpout(maxn) 0209 double precision phxout(maxn) 0210 double precision tmpvec1(maxn), tmpvec2(maxn) 0211 double precision phrwork(3*maxn) a5276edbc9 Patr*0212 cph double precision metricLoc(maxn,maxn) 0213 cph double precision metricTriag(maxn,maxn) 0214 cph double precision metricInv(maxn,maxn) ee5ca1ad15 Patr*0215 0216 cph DATA (metricInv(ii,1),ii=1,maxn) b35bd3101a Jean*0217 cph & / 9.9896D+07, 0.0519D+07, -0.0415D+07, ee5ca1ad15 Patr*0218 cph & 0.0486D+07, -0.2432D+07, 0.1945D+07 / 0219 cph DATA (metricInv(ii,2),ii=1,maxn) 0220 cph & / 0.0519D+07, 9.7403D+07, 0.2077D+07, 0221 cph & -0.2432D+07, 1.2158D+07, -0.9726D+07 / 0222 cph DATA (metricInv(ii,3),ii=1,maxn) 0223 cph & / -0.0415D+07, 0.2077D+07, 9.8338D+07, 0224 cph & 0.1945D+07, -0.9726D+07, 0.7781D+07 / 0225 cph DATA (metricInv(ii,4),ii=1,maxn) 0226 cph & / 0.0486D+07, -0.2432D+07, 0.1945D+07, 0227 cph & 9.7723D+07, 1.1385D+07, -0.9108D+07 / 0228 cph DATA (metricInv(ii,5),ii=1,maxn) 0229 cph & / -0.2432D+07, 1.2158D+07, -0.9726D+07, 0230 cph & 1.1385D+07, 4.3073D+07, 4.5542D+07 / 0231 cph DATA (metricInv(ii,6),ii=1,maxn) 0232 cph & / 0.1945D+07, -0.9726D+07, 0.7781D+07, 0233 cph & -0.9108D+07, 4.5542D+07, 6.3567D+07 / 0234 cph) 0235 c 0236 c %-----------------------% 0237 c | Executable Statements | 0238 c %-----------------------% 0239 c 0240 c %-------------------------------------------------% 0241 c | The following include statement and assignments | 0242 c | initiate trace output from the internal | 0243 c | actions of ARPACK. See debug.doc in the | 0244 c | DOCUMENTS directory for usage. Initially, the | 0245 c | most useful information will be a breakdown of | 0246 c | time spent in the various stages of computation | 0247 c | given by setting msaupd = 1. | 0248 c %-------------------------------------------------% 0249 c 0250 include 'arpack_debug.h' 0251 ndigit = -3 0252 logfil = 6 0253 msgets = 0 b35bd3101a Jean*0254 msaitr = 0 ee5ca1ad15 Patr*0255 msapps = 0 0256 msaupd = 1 0257 msaup2 = 0 0258 mseigt = 0 0259 mseupd = 0 0260 c 0261 c %-------------------------------------------------% 0262 c | The following sets dimensions for this problem. | 0263 c %-------------------------------------------------% 0264 c 0265 cph( e4939baa16 Patr*0266 m = admtlmrec 0267 n = admtlmrec ee5ca1ad15 Patr*0268 cph) 0269 c 0270 c %------------------------------------------------% b35bd3101a Jean*0271 c | Specifications for ARPACK usage are set | ee5ca1ad15 Patr*0272 c | below: | 0273 c | | b35bd3101a Jean*0274 c | 1) NEV = 4 asks for 4 singular values to be | 0275 c | computed. | ee5ca1ad15 Patr*0276 c | | 0277 c | 2) NCV = 20 sets the length of the Arnoldi | 0278 c | factorization | 0279 c | | 0280 c | 3) This is a standard problem | 0281 c | (indicated by bmat = 'I') | 0282 c | | 0283 c | 4) Ask for the NEV singular values of | 0284 c | largest magnitude | 0285 c | (indicated by which = 'LM') | 0286 c | See documentation in DSAUPD for the | b35bd3101a Jean*0287 c | other options SM, BE. | ee5ca1ad15 Patr*0288 c | | 0289 c | Note: NEV and NCV must satisfy the following | 0290 c | conditions: | 0291 c | NEV <= MAXNEV, | 0292 c | NEV + 1 <= NCV <= MAXNCV | 0293 c %------------------------------------------------% 0294 c 0295 cph( 0296 nev = 15 0297 ncv = 30 0298 bmat = 'I' 0299 cph) 0300 which = 'LM' 0301 c 0302 if ( n .gt. maxn ) then 0303 print *, ' ERROR with _SVD: N is greater than MAXN ' 0304 go to 9000 0305 else if ( m .gt. maxm ) then 0306 print *, ' ERROR with _SVD: M is greater than MAXM ' 0307 go to 9000 0308 else if ( nev .gt. maxnev ) then 0309 print *, ' ERROR with _SVD: NEV is greater than MAXNEV ' 0310 go to 9000 0311 else if ( ncv .gt. maxncv ) then 0312 print *, ' ERROR with _SVD: NCV is greater than MAXNCV ' 0313 go to 9000 0314 end if 0315 c 0316 c %-----------------------------------------------------% 0317 c | Specification of stopping rules and initial | 0318 c | conditions before calling DSAUPD | 0319 c | | 0320 c | abs(sigmaC - sigmaT) < TOL*abs(sigmaC) | 0321 c | computed true | 0322 c | | 0323 c | If TOL .le. 0, then TOL <- macheps | 0324 c | (machine precision) is used. | 0325 c | | 0326 c | IDO is the REVERSE COMMUNICATION parameter | 0327 c | used to specify actions to be taken on return | 0328 c | from DSAUPD. (See usage below.) | 0329 c | | 0330 c | It MUST initially be set to 0 before the first | b35bd3101a Jean*0331 c | call to DSAUPD. | ee5ca1ad15 Patr*0332 c | | 0333 c | INFO on entry specifies starting vector information | 0334 c | and on return indicates error codes | 0335 c | | b35bd3101a Jean*0336 c | Initially, setting INFO=0 indicates that a | ee5ca1ad15 Patr*0337 c | random starting vector is requested to | 0338 c | start the ARNOLDI iteration. Setting INFO to | 0339 c | a nonzero value on the initial call is used | 0340 c | if you want to specify your own starting | b35bd3101a Jean*0341 c | vector (This vector must be placed in RESID.) | ee5ca1ad15 Patr*0342 c | | b35bd3101a Jean*0343 c | The work array WORKL is used in DSAUPD as | ee5ca1ad15 Patr*0344 c | workspace. Its dimension LWORKL is set as | 0345 c | illustrated below. | 0346 c %-----------------------------------------------------% 0347 c 0348 lworkl = ncv*(ncv+8) 0349 cph( b35bd3101a Jean*0350 tol = zero ee5ca1ad15 Patr*0351 cph tol = 1.D-10 0352 cph) 0353 info = 0 0354 ido = 0 0355 c 0356 c %---------------------------------------------------% 0357 c | Specification of Algorithm Mode: | 0358 c | | 0359 c | This program uses the exact shift strategy | 0360 c | (indicated by setting IPARAM(1) = 1.) | 0361 c | IPARAM(3) specifies the maximum number of Arnoldi | 0362 c | iterations allowed. Mode 1 of DSAUPD is used | 0363 c | (IPARAM(7) = 1). All these options can be changed | 0364 c | by the user. For details see the documentation in | 0365 c | DSAUPD. | 0366 c %---------------------------------------------------% 0367 c 0368 ishfts = 1 0369 cph( a5276edbc9 Patr*0370 maxitr = 10 ee5ca1ad15 Patr*0371 c maxitr = 5 0372 c 0373 mode = 1 0374 cph) 0375 iparam(1) = ishfts b35bd3101a Jean*0376 c ee5ca1ad15 Patr*0377 iparam(3) = maxitr b35bd3101a Jean*0378 c ee5ca1ad15 Patr*0379 iparam(7) = mode 0380 c 0381 c %------------------------------------------------% 0382 c | M A I N L O O P (Reverse communication loop) | 0383 c %------------------------------------------------% 0384 c 8f0b59c61c Patr*0385 C-- Set model configuration (fixed arrays) 0386 CALL INITIALISE_FIXED( myThid ) ee5ca1ad15 Patr*0387 c 0388 10 continue 0389 c 0390 print *, 'ph----------------------------------------------------' 0391 print *, 'ph----------------------------------------------------' 0392 print *, 'ph----------------------------------------------------' 0393 c 0394 c %---------------------------------------------% b35bd3101a Jean*0395 c | Repeatedly call the routine DSAUPD and take | ee5ca1ad15 Patr*0396 c | actions indicated by parameter IDO until | 0397 c | either convergence is indicated or maxitr | 0398 c | has been exceeded. | 0399 c %---------------------------------------------% 0400 c faf44775ba Patr*0401 ctest( b35bd3101a Jean*0402 CALL DSAUPD ( ido, bmat, n, which, nev, tol, resid, ee5ca1ad15 Patr*0403 & ncv, v, ldv, iparam, ipntr, workd, workl, 0404 & lworkl, info ) faf44775ba Patr*0405 ctest ido = -1 0406 ctest) ee5ca1ad15 Patr*0407 c 0408 cph( b35bd3101a Jean*0409 print *, 'ph-count: optimcycle, ido, info, ipntr(1) ', a5276edbc9 Patr*0410 & optimcycle, ido, info, ipntr(1) ee5ca1ad15 Patr*0411 cph) 0412 0413 if (ido .eq. -1 .or. ido .eq. 1) then 0414 c 0415 c %---------------------------------------% 0416 c | Perform matrix vector multiplications | 0417 c | w <--- A*x (av()) | 0418 c | y <--- A'*w (atv()) | 0419 c | The user should supply his/her own | 0420 c | matrix vector multiplication routines | 0421 c | here that takes workd(ipntr(1)) as | 0422 c | the input, and returns the result in | 0423 c | workd(ipntr(2)). | 0424 c %---------------------------------------% 0425 c b35bd3101a Jean*0426 c call av (m, n, workd(ipntr(1)), ax) a5276edbc9 Patr*0427 c call atv (m, n, ax, workd(ipntr(2))) ee5ca1ad15 Patr*0428 c 0429 cph( a5276edbc9 Patr*0430 c do l = 1, n b35bd3101a Jean*0431 c print *, 'ph-test A ', l, a5276edbc9 Patr*0432 c & workd(ipntr(1)+l), workd(ipntr(2)+l), workd(l) 0433 c enddo ee5ca1ad15 Patr*0434 cph) 0435 do l = 1, n 875eeaf325 Patr*0436 phtmpadmtlm(l) = workd(ipntr(1)+l-1) ee5ca1ad15 Patr*0437 enddo a5276edbc9 Patr*0438 faf44775ba Patr*0439 IF (optimcycle .GT. 0 ) THEN 0440 _BEGIN_MASTER( mythid ) 0441 IF ( myProcId .eq. 0 ) THEN 0442 call admtlm_dsvd2model( .FALSE., mythid ) 0443 ENDIF 0444 _END_MASTER( mythid ) 0445 cph CALL ADMTLM_UPXX( mythid ) 0446 ENDIF a5276edbc9 Patr*0447 faf44775ba Patr*0448 c-- MWMWMWMWMWMWMW 8f0b59c61c Patr*0449 CALL ADMTLM_DRIVER( mythid ) faf44775ba Patr*0450 c-- MWMWMWMWMWMWMW ee5ca1ad15 Patr*0451 faf44775ba Patr*0452 _BEGIN_MASTER( mythid ) 0453 IF ( myProcId .eq. 0 ) THEN 0454 call admtlm_model2dsvd( .FALSE., mythid ) 0455 ENDIF 0456 _END_MASTER( mythid ) a5276edbc9 Patr*0457 ee5ca1ad15 Patr*0458 do l = 1, n 875eeaf325 Patr*0459 workd(ipntr(2)+l-1) = phtmpadmtlm(l) ee5ca1ad15 Patr*0460 enddo 0461 c a5276edbc9 Patr*0462 faf44775ba Patr*0463 if (optimcycle .EQ. 4) a5276edbc9 Patr*0464 & STOP 'in ADMTLM_DSVD after the_model_main' 0465 ee5ca1ad15 Patr*0466 cph( 0467 cph Since we solve for a generalized EV problem, we have to solve 0468 cph M*y=A*x for y with known matrix M and vector A*x a5276edbc9 Patr*0469 c do ii = 1, n 0470 c phxout(ii) = 0. 0471 c do jj = 1, n 0472 c phxout(ii) = phxout(ii) + 0473 c & metricInv(ii,jj)*phtmpout(jj) 0474 c enddo 0475 c print *, 'ph-test C ', ii, phxout(ii) 0476 c enddo 0477 c ee5ca1ad15 Patr*0478 c call dposvx ( 0479 c & 'E', 'U', maxn, 1, metricLoc, 6, metricTriag, 6, 0480 c & 'N', tmpvec1, phtmpout, 6, phxout, 6, tmpvec2, 0481 c & ferr, berr, phrwork, phiwork, ilinsysinfo) 0482 cph b35bd3101a Jean*0483 c call dgesv ( maxn, 1, metricLoc, maxn, ee5ca1ad15 Patr*0484 c & ipiv, phtmpout, maxn, ilinsysinfo ) 0485 cph 0486 c 0487 cph Finally schift result y -> workd(ipntr(2)) 0488 cph call dcopy ( maxn, tmpvec1, 1, workd(ipntr(2)), 1 ) 0489 cph We have restored orig. standard EV problem 0490 cph OP*x = lambda*x for OP = INV(M)*A 0491 cph) 0492 c 0493 c %-----------------------------------------% 0494 c | L O O P B A C K to call DSAUPD again. | 0495 c %-----------------------------------------% 0496 c a5276edbc9 Patr*0497 optimcycle = optimcycle + 1 0498 ee5ca1ad15 Patr*0499 go to 10 0500 c b35bd3101a Jean*0501 end if ee5ca1ad15 Patr*0502 c 0503 cph( 0504 1001 continue 0505 print *, 'ph-continue ', info, ierr 0506 cph) 0507 c 0508 c %----------------------------------------% 0509 c | Either we have convergence or there is | 0510 c | an error. | 0511 c %----------------------------------------% 0512 c 0513 if ( info .lt. 0 ) then 0514 c 0515 c %--------------------------% 0516 c | Error message. Check the | 0517 c | documentation in DSAUPD. | 0518 c %--------------------------% 0519 c 0520 print *, ' ' 0521 print *, ' Error with _saupd, info = ', info 0522 print *, ' Check documentation in _saupd ' 0523 print *, ' ' 0524 c b35bd3101a Jean*0525 else ee5ca1ad15 Patr*0526 c 0527 c %--------------------------------------------% 0528 c | No fatal errors occurred. | 0529 c | Post-Process using DSEUPD. | 0530 c | | b35bd3101a Jean*0531 c | Computed singular values may be extracted. | ee5ca1ad15 Patr*0532 c | | 0533 c | Singular vectors may also be computed now | b35bd3101a Jean*0534 c | if desired. (indicated by rvec = .true.) | ee5ca1ad15 Patr*0535 c | | 0536 c | The routine DSEUPD now called to do this | b35bd3101a Jean*0537 c | post processing | ee5ca1ad15 Patr*0538 c %--------------------------------------------% b35bd3101a Jean*0539 c ee5ca1ad15 Patr*0540 rvec = .true. 0541 c b35bd3101a Jean*0542 call dseupd ( rvec, 'All', select, s, v, ldv, sigma, 0543 & bmat, n, which, nev, tol, resid, ncv, v, ldv, ee5ca1ad15 Patr*0544 & iparam, ipntr, workd, workl, lworkl, ierr ) 0545 c 0546 c %-----------------------------------------------% 0547 c | Singular values are returned in the first | 0548 c | column of the two dimensional array S | b35bd3101a Jean*0549 c | and the corresponding right singular vectors | ee5ca1ad15 Patr*0550 c | are returned in the first NEV columns of the | 0551 c | two dimensional array V as requested here. | 0552 c %-----------------------------------------------% 0553 c 0554 if ( ierr .ne. 0) then 0555 c 0556 c %------------------------------------% 0557 c | Error condition: | 0558 c | Check the documentation of DSEUPD. | 0559 c %------------------------------------% 0560 c 0561 print *, ' ' 0562 print *, ' Error with _seupd, info = ', ierr 0563 print *, ' Check the documentation of _seupd. ' 0564 print *, ' ' 0565 c 0566 else 0567 c 0568 nconv = iparam(5) 0569 cph( 0570 do j=1, nconv 0571 print *, 'ph-ev ', j, s(j,1), s(j,2) 0572 enddo 0573 0574 STOP 'TEST AFTER dneupd' 0575 cph) 0576 do 20 j=1, nconv 0577 c 0578 s(j,1) = sqrt(s(j,1)) 0579 c 0580 c %-----------------------------% 0581 c | Compute the left singular | 0582 c | vectors from the formula | 0583 c | | 0584 c | u = Av/sigma | 0585 c | | 0586 c | u should have norm 1 so | 0587 c | divide by norm(Av) instead. | 0588 c %-----------------------------% 0589 c 0590 do l = 1, n 0591 tmpvec1(l) = v(l,j) 0592 enddo a5276edbc9 Patr*0593 ee5ca1ad15 Patr*0594 c 0595 c-- MWMWMWMWMWMWMW a5276edbc9 Patr*0596 cph call box_main( tmpvec1, tmpvec2, metricLoc, ldoadjoint ) ee5ca1ad15 Patr*0597 c-- MWMWMWMWMWMWMW 0598 c 0599 call dcopy( m, tmpvec2, 1, u(l,j), 1 ) 0600 temp = one / dnrm2( m, u(l,j), 1 ) 0601 cph( 0602 print *, 'ph-print ', j, dnrm2( m, u(l,j), 1 ) 0603 cph) 0604 call dscal(m , temp, u(l,j), 1 ) 0605 0606 cph do l = 1, n 0607 cph u(l,j) = tmpvec2(l) 0608 cph enddo 0609 cph temp = one/dnrm2(m, tmpvec2, 1) 0610 cph call dscal(m, temp, tmpvec2, 1) 0611 c 0612 c %---------------------------% 0613 c | | 0614 c | Compute the residual norm | 0615 c | | 0616 c | || A*v - sigma*u || | 0617 c | | 0618 c | for the NCONV accurately | 0619 c | computed singular values | 0620 c | and vectors. (iparam(5) | 0621 c | indicates how many are | 0622 c | accurate to the requested | 0623 c | tolerance). | 0624 c | Store the result in 2nd | 0625 c | column of array S. | 0626 c %---------------------------% 0627 c 0628 call daxpy(m, -s(j,1), u(1,j), 1, tmpvec2, 1) 0629 s(j,2) = dnrm2(m, tmpvec2, 1) 0630 c 0631 20 continue 0632 c 0633 c %-------------------------------% 0634 c | Display computed residuals | 0635 c %-------------------------------% 0636 c 0637 call dmout(6, nconv, 2, s, maxncv, -6, 0638 & 'Singular values and direct residuals') 0639 end if 0640 c 0641 c %------------------------------------------% 0642 c | Print additional convergence information | 0643 c %------------------------------------------% 0644 c 0645 if ( info .eq. 1) then 0646 print *, ' ' 0647 print *, ' Maximum number of iterations reached.' 0648 print *, ' ' 0649 else if ( info .eq. 3) then b35bd3101a Jean*0650 print *, ' ' ee5ca1ad15 Patr*0651 print *, ' No shifts could be applied during implicit', 0652 & ' Arnoldi update, try increasing NCV.' 0653 print *, ' ' b35bd3101a Jean*0654 end if ee5ca1ad15 Patr*0655 c 0656 print *, ' ' 0657 print *, ' _SVD ' 0658 print *, ' ==== ' 0659 print *, ' ' 0660 print *, ' Size of the matrix is ', n 0661 print *, ' The number of Ritz values requested is ', nev 0662 print *, ' The number of Arnoldi vectors generated', 0663 & ' (NCV) is ', ncv 0664 print *, ' What portion of the spectrum: ', which b35bd3101a Jean*0665 print *, ' The number of converged Ritz values is ', 0666 & nconv ee5ca1ad15 Patr*0667 print *, ' The number of Implicit Arnoldi update', 0668 & ' iterations taken is ', iparam(3) 0669 print *, ' The number of OP*x is ', iparam(9) 0670 print *, ' The convergence criterion is ', tol 0671 print *, ' ' 0672 c 0673 end if 0674 c 0675 c %-------------------------% 0676 c | Done with program dsvd. | 0677 c %-------------------------% 0678 c 0679 9000 continue 0680 c 0681 end
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine. The LXR team |
|