Branch data Line data Source code
1 : : !!****m* ABINIT/m_hybrd
2 : : !! NAME
3 : : !! m_hybrd
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains a modified Powell method for root-finding.
7 : : !! This was taken from https://www.netlib.org/minpack/.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (c) 2012 Ondrej Certik
11 : : !!
12 : : !! Minpack Copyright Notice (1999) University of Chicago. All rights reserved
13 : : !!
14 : : !! Redistribution and use in source and binary forms, with or
15 : : !! without modification, are permitted provided that the
16 : : !! following conditions are met:
17 : : !!
18 : : !! 1. Redistributions of source code must retain the above
19 : : !! copyright notice, this list of conditions and the following
20 : : !! disclaimer.
21 : : !!
22 : : !! 2. Redistributions in binary form must reproduce the above
23 : : !! copyright notice, this list of conditions and the following
24 : : !! disclaimer in the documentation and/or other materials
25 : : !! provided with the distribution.
26 : : !!
27 : : !! 3. The end-user documentation included with the
28 : : !! redistribution, if any, must include the following
29 : : !! acknowledgment:
30 : : !!
31 : : !! "This product includes software developed by the
32 : : !! University of Chicago, as Operator of Argonne National
33 : : !! Laboratory.
34 : : !!
35 : : !! Alternately, this acknowledgment may appear in the software
36 : : !! itself, if and wherever such third-party acknowledgments
37 : : !! normally appear.
38 : : !!
39 : : !! 4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
40 : : !! WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
41 : : !! UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
42 : : !! THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
43 : : !! IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
44 : : !! OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
45 : : !! OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
46 : : !! OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
47 : : !! USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
48 : : !! THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
49 : : !! DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
50 : : !! UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
51 : : !! BE CORRECTED.
52 : : !!
53 : : !! 5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
54 : : !! HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
55 : : !! ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
56 : : !! INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
57 : : !! ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
58 : : !! PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
59 : : !! SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
60 : : !! (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
61 : : !! EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
62 : : !! POSSIBILITY OF SUCH LOSS OR DAMAGES.
63 : : !!
64 : : !! This file is distributed under the terms of the
65 : : !! GNU General Public License, see ~abinit/COPYING
66 : : !! or http://www.gnu.org/copyleft/gpl.txt .
67 : : !!
68 : : !! SOURCE
69 : :
70 : : #if defined HAVE_CONFIG_H
71 : : #include "config.h"
72 : : #endif
73 : :
74 : : MODULE m_hybrd
75 : :
76 : : implicit none
77 : :
78 : : private
79 : :
80 : : public :: hybrd
81 : :
82 : : CONTAINS !========================================================================================
83 : : !!***
84 : :
85 : : recursive &
86 : 0 : & subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag, &
87 : 0 : & mode,factor,nprint,info,nfev,fjac,ldfjac,r,lr, &
88 : 0 : & qtf,wa1,wa2,wa3,wa4)
89 : : integer n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,lr
90 : : double precision xtol,epsfcn,factor
91 : : double precision x(n),fvec(n),diag(n),fjac(ldfjac,n),r(lr), &
92 : : & qtf(n),wa1(n),wa2(n),wa3(n),wa4(n)
93 : : interface
94 : : subroutine fcn(n,x,fvec,iflag)
95 : : integer, intent(in) :: n,iflag
96 : : double precision, intent(in) :: x(n)
97 : : double precision, intent(inout) :: fvec(n)
98 : : end subroutine fcn
99 : : end interface
100 : : !c **********
101 : : !c
102 : : !c subroutine hybrd
103 : : !c
104 : : !c the purpose of hybrd is to find a zero of a system of
105 : : !c n nonlinear functions in n variables by a modification
106 : : !c of the powell hybrid method. the user must provide a
107 : : !c subroutine which calculates the functions. the jacobian is
108 : : !c then calculated by a forward-difference approximation.
109 : : !c
110 : : !c the subroutine statement is
111 : : !c
112 : : !c subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,
113 : : !c diag,mode,factor,nprint,info,nfev,fjac,
114 : : !c ldfjac,r,lr,qtf,wa1,wa2,wa3,wa4)
115 : : !c
116 : : !c where
117 : : !c
118 : : !c fcn is the name of the user-supplied subroutine which
119 : : !c calculates the functions. fcn must be declared
120 : : !c in an external statement in the user calling
121 : : !c program, and should be written as follows.
122 : : !c
123 : : !c subroutine fcn(n,x,fvec,iflag)
124 : : !c integer n,iflag
125 : : !c double precision x(n),fvec(n)
126 : : !c ----------
127 : : !c calculate the functions at x and
128 : : !c return this vector in fvec.
129 : : !c ---------
130 : : !c return
131 : : !c end
132 : : !c
133 : : !c the value of iflag should not be changed by fcn unless
134 : : !c the user wants to terminate execution of hybrd.
135 : : !c in this case set iflag to a negative integer.
136 : : !c
137 : : !c n is a positive integer input variable set to the number
138 : : !c of functions and variables.
139 : : !c
140 : : !c x is an array of length n. on input x must contain
141 : : !c an initial estimate of the solution vector. on output x
142 : : !c contains the final estimate of the solution vector.
143 : : !c
144 : : !c fvec is an output array of length n which contains
145 : : !c the functions evaluated at the output x.
146 : : !c
147 : : !c xtol is a nonnegative input variable. termination
148 : : !c occurs when the relative error between two consecutive
149 : : !c iterates is at most xtol.
150 : : !c
151 : : !c maxfev is a positive integer input variable. termination
152 : : !c occurs when the number of calls to fcn is at least maxfev
153 : : !c by the end of an iteration.
154 : : !c
155 : : !c ml is a nonnegative integer input variable which specifies
156 : : !c the number of subdiagonals within the band of the
157 : : !c jacobian matrix. if the jacobian is not banded, set
158 : : !c ml to at least n - 1.
159 : : !c
160 : : !c mu is a nonnegative integer input variable which specifies
161 : : !c the number of superdiagonals within the band of the
162 : : !c jacobian matrix. if the jacobian is not banded, set
163 : : !c mu to at least n - 1.
164 : : !c
165 : : !c epsfcn is an input variable used in determining a suitable
166 : : !c step length for the forward-difference approximation. this
167 : : !c approximation assumes that the relative errors in the
168 : : !c functions are of the order of epsfcn. if epsfcn is less
169 : : !c than the machine precision, it is assumed that the relative
170 : : !c errors in the functions are of the order of the machine
171 : : !c precision.
172 : : !c
173 : : !c diag is an array of length n. if mode = 1 (see
174 : : !c below), diag is internally set. if mode = 2, diag
175 : : !c must contain positive entries that serve as
176 : : !c multiplicative scale factors for the variables.
177 : : !c
178 : : !c mode is an integer input variable. if mode = 1, the
179 : : !c variables will be scaled internally. if mode = 2,
180 : : !c the scaling is specified by the input diag. other
181 : : !c values of mode are equivalent to mode = 1.
182 : : !c
183 : : !c factor is a positive input variable used in determining the
184 : : !c initial step bound. this bound is set to the product of
185 : : !c factor and the euclidean norm of diag*x if nonzero, or else
186 : : !c to factor itself. in most cases factor should lie in the
187 : : !c interval (.1,100.). 100. is a generally recommended value.
188 : : !c
189 : : !c nprint is an integer input variable that enables controlled
190 : : !c printing of iterates if it is positive. in this case,
191 : : !c fcn is called with iflag = 0 at the beginning of the first
192 : : !c iteration and every nprint iterations thereafter and
193 : : !c immediately prior to return, with x and fvec available
194 : : !c for printing. if nprint is not positive, no special calls
195 : : !c of fcn with iflag = 0 are made.
196 : : !c
197 : : !c info is an integer output variable. if the user has
198 : : !c terminated execution, info is set to the (negative)
199 : : !c value of iflag. see description of fcn. otherwise,
200 : : !c info is set as follows.
201 : : !c
202 : : !c info = 0 improper input parameters.
203 : : !c
204 : : !c info = 1 relative error between two consecutive iterates
205 : : !c is at most xtol.
206 : : !c
207 : : !c info = 2 number of calls to fcn has reached or exceeded
208 : : !c maxfev.
209 : : !c
210 : : !c info = 3 xtol is too small. no further improvement in
211 : : !c the approximate solution x is possible.
212 : : !c
213 : : !c info = 4 iteration is not making good progress, as
214 : : !c measured by the improvement from the last
215 : : !c five jacobian evaluations.
216 : : !c
217 : : !c info = 5 iteration is not making good progress, as
218 : : !c measured by the improvement from the last
219 : : !c ten iterations.
220 : : !c
221 : : !c nfev is an integer output variable set to the number of
222 : : !c calls to fcn.
223 : : !c
224 : : !c fjac is an output n by n array which contains the
225 : : !c orthogonal matrix q produced by the qr factorization
226 : : !c of the final approximate jacobian.
227 : : !c
228 : : !c ldfjac is a positive integer input variable not less than n
229 : : !c which specifies the leading dimension of the array fjac.
230 : : !c
231 : : !c r is an output array of length lr which contains the
232 : : !c upper triangular matrix produced by the qr factorization
233 : : !c of the final approximate jacobian, stored rowwise.
234 : : !c
235 : : !c lr is a positive integer input variable not less than
236 : : !c (n*(n+1))/2.
237 : : !c
238 : : !c qtf is an output array of length n which contains
239 : : !c the vector (q transpose)*fvec.
240 : : !c
241 : : !c wa1, wa2, wa3, and wa4 are work arrays of length n.
242 : : !c
243 : : !c subprograms called
244 : : !c
245 : : !c user-supplied ...... fcn
246 : : !c
247 : : !c minpack-supplied ... dogleg,dpmpar,enorm,fdjac1,
248 : : !c qform,qrfac,r1mpyq,r1updt
249 : : !c
250 : : !c fortran-supplied ... dabs,dmax1,dmin1,min0,mod
251 : : !c
252 : : !c argonne national laboratory. minpack project. march 1980.
253 : : !c burton s. garbow, kenneth e. hillstrom, jorge j. more
254 : : !c
255 : : !c **********
256 : : integer i,iflag,iter,j,jm1,l,msum,ncfail,ncsuc,nslow1,nslow2
257 : : integer iwa(1)
258 : : logical jeval,sing
259 : : double precision actred,delta,epsmch,fnorm,fnorm1,one,pnorm, &
260 : : & prered,p1,p5,p001,p0001,ratio,sum,temp,xnorm, &
261 : : & zero
262 : : !double precision dpmpar,enorm
263 : : data one,p1,p5,p001,p0001,zero &
264 : : & /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/
265 : : !c
266 : : !c epsmch is the machine precision.
267 : : !c
268 : 0 : epsmch = dpmpar(1)
269 : : !c
270 : 0 : info = 0
271 : 0 : iflag = 0
272 : 0 : nfev = 0
273 : : !c
274 : : !c check the input parameters for errors.
275 : : !c
276 : : if (n .le. 0 .or. xtol .lt. zero .or. maxfev .le. 0 &
277 : : & .or. ml .lt. 0 .or. mu .lt. 0 .or. factor .le. zero &
278 [ # # # # : 0 : & .or. ldfjac .lt. n .or. lr .lt. (n*(n + 1))/2) go to 300
# # # # #
# # # #
# ]
279 [ # # ]: 0 : if (mode .ne. 2) go to 20
280 [ # # ]: 0 : do 10 j = 1, n
281 [ # # ]: 0 : if (diag(j) .le. zero) go to 300
282 : 0 : 10 continue
283 : : 20 continue
284 : : !c
285 : : !c evaluate the function at the starting point
286 : : !c and calculate its norm.
287 : : !c
288 : 0 : iflag = 1
289 : 0 : call fcn(n,x,fvec,iflag)
290 : 0 : nfev = 1
291 : : if (iflag .lt. 0) go to 300
292 : 0 : fnorm = enorm(n,fvec)
293 : : !c
294 : : !c determine the number of calls to fcn needed to compute
295 : : !c the jacobian matrix.
296 : : !c
297 : 0 : msum = min0(ml+mu+1,n)
298 : : !c
299 : : !c initialize iteration counter and monitors.
300 : : !c
301 : 0 : iter = 1
302 : 0 : ncsuc = 0
303 : 0 : ncfail = 0
304 : 0 : nslow1 = 0
305 : 0 : nslow2 = 0
306 : : !c
307 : : !c beginning of the outer loop.
308 : : !c
309 : : 30 continue
310 : 0 : jeval = .true.
311 : : !c
312 : : !c calculate the jacobian matrix.
313 : : !c
314 : 0 : iflag = 2
315 : : call fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1, &
316 : 0 : & wa2)
317 : 0 : nfev = nfev + msum
318 : : if (iflag .lt. 0) go to 300
319 : : !c
320 : : !c compute the qr factorization of the jacobian.
321 : : !c
322 : 0 : call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
323 : : !c
324 : : !c on the first iteration and if mode is 1, scale according
325 : : !c to the norms of the columns of the initial jacobian.
326 : : !c
327 [ # # ]: 0 : if (iter .ne. 1) go to 70
328 [ # # ]: 0 : if (mode .eq. 2) go to 50
329 [ # # ]: 0 : do 40 j = 1, n
330 : 0 : diag(j) = wa2(j)
331 [ # # ]: 0 : if (wa2(j) .eq. zero) diag(j) = one
332 : 0 : 40 continue
333 : : 50 continue
334 : : !c
335 : : !c on the first iteration, calculate the norm of the scaled x
336 : : !c and initialize the step bound delta.
337 : : !c
338 [ # # ]: 0 : do 60 j = 1, n
339 : 0 : wa3(j) = diag(j)*x(j)
340 : 0 : 60 continue
341 : 0 : xnorm = enorm(n,wa3)
342 : 0 : delta = factor*xnorm
343 [ # # ]: 0 : if (delta .eq. zero) delta = factor
344 : : 70 continue
345 : : !c
346 : : !c form (q transpose)*fvec and store in qtf.
347 : : !c
348 [ # # ]: 0 : do 80 i = 1, n
349 : 0 : qtf(i) = fvec(i)
350 : 0 : 80 continue
351 [ # # ]: 0 : do 120 j = 1, n
352 [ # # ]: 0 : if (fjac(j,j) .eq. zero) go to 110
353 : : sum = zero
354 [ # # ]: 0 : do 90 i = j, n
355 : 0 : sum = sum + fjac(i,j)*qtf(i)
356 : 0 : 90 continue
357 : 0 : temp = -sum/fjac(j,j)
358 [ # # ]: 0 : do 100 i = j, n
359 : 0 : qtf(i) = qtf(i) + fjac(i,j)*temp
360 : 0 : 100 continue
361 : : 110 continue
362 : 0 : 120 continue
363 : : !c
364 : : !c copy the triangular factor of the qr factorization into r.
365 : : !c
366 : 0 : sing = .false.
367 [ # # ]: 0 : do 150 j = 1, n
368 : 0 : l = j
369 : 0 : jm1 = j - 1
370 [ # # ]: 0 : if (jm1 .lt. 1) go to 140
371 [ # # ]: 0 : do 130 i = 1, jm1
372 : 0 : r(l) = fjac(i,j)
373 : 0 : l = l + n - i
374 : 0 : 130 continue
375 : : 140 continue
376 : 0 : r(l) = wa1(j)
377 [ # # ]: 0 : if (wa1(j) .eq. zero) sing = .true.
378 : 0 : 150 continue
379 : : !c
380 : : !c accumulate the orthogonal factor in fjac.
381 : : !c
382 : 0 : call qform(n,n,fjac,ldfjac,wa1)
383 : : !c
384 : : !c rescale if necessary.
385 : : !c
386 [ # # ]: 0 : if (mode .eq. 2) go to 170
387 [ # # ]: 0 : do 160 j = 1, n
388 : 0 : diag(j) = dmax1(diag(j),wa2(j))
389 : 0 : 160 continue
390 : : 170 continue
391 : : !c
392 : : !c beginning of the inner loop.
393 : : !c
394 : : 180 continue
395 : : !c
396 : : !c if requested, call fcn to enable printing of iterates.
397 : : !c
398 [ # # ]: 0 : if (nprint .le. 0) go to 190
399 : 0 : iflag = 0
400 [ # # ]: 0 : if (mod(iter-1,nprint) .eq. 0) call fcn(n,x,fvec,iflag)
401 : : if (iflag .lt. 0) go to 300
402 : : 190 continue
403 : : !c
404 : : !c determine the direction p.
405 : : !c
406 : 0 : call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
407 : : !c
408 : : !c store the direction p and x + p. calculate the norm of p.
409 : : !c
410 [ # # ]: 0 : do 200 j = 1, n
411 : 0 : wa1(j) = -wa1(j)
412 : 0 : wa2(j) = x(j) + wa1(j)
413 : 0 : wa3(j) = diag(j)*wa1(j)
414 : 0 : 200 continue
415 : 0 : pnorm = enorm(n,wa3)
416 : : !c
417 : : !c on the first iteration, adjust the initial step bound.
418 : : !c
419 [ # # ]: 0 : if (iter .eq. 1) delta = dmin1(delta,pnorm)
420 : : !c
421 : : !c evaluate the function at x + p and calculate its norm.
422 : : !c
423 : 0 : iflag = 1
424 : 0 : call fcn(n,wa2,wa4,iflag)
425 : 0 : nfev = nfev + 1
426 : : if (iflag .lt. 0) go to 300
427 : 0 : fnorm1 = enorm(n,wa4)
428 : : !c
429 : : !c compute the scaled actual reduction.
430 : : !c
431 : 0 : actred = -one
432 [ # # ]: 0 : if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
433 : : !c
434 : : !c compute the scaled predicted reduction.
435 : : !c
436 : 0 : l = 1
437 [ # # ]: 0 : do 220 i = 1, n
438 : 0 : sum = zero
439 [ # # ]: 0 : do 210 j = i, n
440 : 0 : sum = sum + r(l)*wa1(j)
441 : 0 : l = l + 1
442 : 0 : 210 continue
443 : 0 : wa3(i) = qtf(i) + sum
444 : 0 : 220 continue
445 : 0 : temp = enorm(n,wa3)
446 : 0 : prered = zero
447 [ # # ]: 0 : if (temp .lt. fnorm) prered = one - (temp/fnorm)**2
448 : : !c
449 : : !c compute the ratio of the actual to the predicted
450 : : !c reduction.
451 : : !c
452 : 0 : ratio = zero
453 [ # # ]: 0 : if (prered .gt. zero) ratio = actred/prered
454 : : !c
455 : : !c update the step bound.
456 : : !c
457 [ # # ]: 0 : if (ratio .ge. p1) go to 230
458 : 0 : ncsuc = 0
459 : 0 : ncfail = ncfail + 1
460 : 0 : delta = p5*delta
461 : 0 : go to 240
462 : : 230 continue
463 : 0 : ncfail = 0
464 : 0 : ncsuc = ncsuc + 1
465 [ # # # # ]: 0 : if (ratio .ge. p5 .or. ncsuc .gt. 1) &
466 : 0 : & delta = dmax1(delta,pnorm/p5)
467 [ # # ]: 0 : if (dabs(ratio-one) .le. p1) delta = pnorm/p5
468 : : 240 continue
469 : : !c
470 : : !c test for successful iteration.
471 : : !c
472 [ # # ]: 0 : if (ratio .lt. p0001) go to 260
473 : : !c
474 : : !c successful iteration. update x, fvec, and their norms.
475 : : !c
476 [ # # ]: 0 : do 250 j = 1, n
477 : 0 : x(j) = wa2(j)
478 : 0 : wa2(j) = diag(j)*x(j)
479 : 0 : fvec(j) = wa4(j)
480 : 0 : 250 continue
481 : 0 : xnorm = enorm(n,wa2)
482 : 0 : fnorm = fnorm1
483 : 0 : iter = iter + 1
484 : : 260 continue
485 : : !c
486 : : !c determine the progress of the iteration.
487 : : !c
488 : 0 : nslow1 = nslow1 + 1
489 [ # # ]: 0 : if (actred .ge. p001) nslow1 = 0
490 [ # # ]: 0 : if (jeval) nslow2 = nslow2 + 1
491 [ # # ]: 0 : if (actred .ge. p1) nslow2 = 0
492 : : !c
493 : : !c test for convergence.
494 : : !c
495 [ # # # # ]: 0 : if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1
496 [ # # ]: 0 : if (info .ne. 0) go to 300
497 : : !c
498 : : !c tests for termination and stringent tolerances.
499 : : !c
500 [ # # ]: 0 : if (nfev .ge. maxfev) info = 2
501 [ # # ]: 0 : if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3
502 [ # # ]: 0 : if (nslow2 .eq. 5) info = 4
503 [ # # ]: 0 : if (nslow1 .eq. 10) info = 5
504 [ # # ]: 0 : if (info .ne. 0) go to 300
505 : : !c
506 : : !c criterion for recalculating jacobian approximation
507 : : !c by forward differences.
508 : : !c
509 [ # # ]: 0 : if (ncfail .eq. 2) go to 290
510 : : !c
511 : : !c calculate the rank one modification to the jacobian
512 : : !c and update qtf if necessary.
513 : : !c
514 [ # # ]: 0 : do 280 j = 1, n
515 : : sum = zero
516 [ # # ]: 0 : do 270 i = 1, n
517 : 0 : sum = sum + fjac(i,j)*wa4(i)
518 : 0 : 270 continue
519 : 0 : wa2(j) = (sum - wa3(j))/pnorm
520 : 0 : wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
521 [ # # ]: 0 : if (ratio .ge. p0001) qtf(j) = sum
522 : 0 : 280 continue
523 : : !c
524 : : !c compute the qr factorization of the updated jacobian.
525 : : !c
526 : 0 : call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
527 : 0 : call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
528 : 0 : call r1mpyq(1,n,qtf,1,wa2,wa3)
529 : : !c
530 : : !c end of the inner loop.
531 : : !c
532 : 0 : jeval = .false.
533 : 0 : go to 180
534 : : 290 continue
535 : : !c
536 : : !c end of the outer loop.
537 : : !c
538 : 0 : go to 30
539 : : 300 continue
540 : : !c
541 : : !c termination, either normal or user imposed.
542 : : !c
543 [ # # ]: 0 : if (iflag .lt. 0) info = iflag
544 : 0 : iflag = 0
545 [ # # ]: 0 : if (nprint .gt. 0) call fcn(n,x,fvec,iflag)
546 : 0 : return
547 : : !c
548 : : !c last card of subroutine hybrd.
549 : : !c
550 : : end subroutine hybrd
551 : :
552 : : recursive &
553 : 0 : &double precision function dpmpar(i)
554 : : integer i
555 : : !c **********
556 : : !c
557 : : !c Function dpmpar
558 : : !c
559 : : !c This function provides double precision machine parameters
560 : : !c when the appropriate set of data statements is activated (by
561 : : !c removing the c from column 1) and all other data statements are
562 : : !c rendered inactive. Most of the parameter values were obtained
563 : : !c from the corresponding Bell Laboratories Port Library function.
564 : : !c
565 : : !c The function statement is
566 : : !c
567 : : !c double precision function dpmpar(i)
568 : : !c
569 : : !c where
570 : : !c
571 : : !c i is an integer input variable set to 1, 2, or 3 which
572 : : !c selects the desired machine parameter. If the machine has
573 : : !c t base b digits and its smallest and largest exponents are
574 : : !c emin and emax, respectively, then these parameters are
575 : : !c
576 : : !c dpmpar(1) = b**(1 - t), the machine precision,
577 : : !c
578 : : !c dpmpar(2) = b**(emin - 1), the smallest magnitude,
579 : : !c
580 : : !c dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude.
581 : : !c
582 : : !c Argonne National Laboratory. MINPACK Project. November 1996.
583 : : !c Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More'
584 : : !c
585 : : !c **********
586 : : integer mcheps(4)
587 : : integer minmag(4)
588 : : integer maxmag(4)
589 : : double precision dmach(3)
590 : : equivalence (dmach(1),mcheps(1))
591 : : equivalence (dmach(2),minmag(1))
592 : : equivalence (dmach(3),maxmag(1))
593 : : !c
594 : : !c Machine constants for the IBM 360/370 series,
595 : : !c the Amdahl 470/V6, the ICL 2900, the Itel AS/6,
596 : : !c the Xerox Sigma 5/7/9 and the Sel systems 85/86.
597 : : !c
598 : : !c data mcheps(1),mcheps(2) / z34100000, z00000000 /
599 : : !c data minmag(1),minmag(2) / z00100000, z00000000 /
600 : : !c data maxmag(1),maxmag(2) / z7fffffff, zffffffff /
601 : : !c
602 : : !c Machine constants for the Honeywell 600/6000 series.
603 : : !c
604 : : !c data mcheps(1),mcheps(2) / o606400000000, o000000000000 /
605 : : !c data minmag(1),minmag(2) / o402400000000, o000000000000 /
606 : : !c data maxmag(1),maxmag(2) / o376777777777, o777777777777 /
607 : : !c
608 : : !c Machine constants for the CDC 6000/7000 series.
609 : : !c
610 : : !c data mcheps(1) / 15614000000000000000b /
611 : : !c data mcheps(2) / 15010000000000000000b /
612 : : !c
613 : : !c data minmag(1) / 00604000000000000000b /
614 : : !c data minmag(2) / 00000000000000000000b /
615 : : !c
616 : : !c data maxmag(1) / 37767777777777777777b /
617 : : !c data maxmag(2) / 37167777777777777777b /
618 : : !c
619 : : !c Machine constants for the PDP-10 (KA processor).
620 : : !c
621 : : !c data mcheps(1),mcheps(2) / "114400000000, "000000000000 /
622 : : !c data minmag(1),minmag(2) / "033400000000, "000000000000 /
623 : : !c data maxmag(1),maxmag(2) / "377777777777, "344777777777 /
624 : : !c
625 : : !c Machine constants for the PDP-10 (KI processor).
626 : : !c
627 : : !c data mcheps(1),mcheps(2) / "104400000000, "000000000000 /
628 : : !c data minmag(1),minmag(2) / "000400000000, "000000000000 /
629 : : !c data maxmag(1),maxmag(2) / "377777777777, "377777777777 /
630 : : !c
631 : : !c Machine constants for the PDP-11.
632 : : !c
633 : : !c data mcheps(1),mcheps(2) / 9472, 0 /
634 : : !c data mcheps(3),mcheps(4) / 0, 0 /
635 : : !c
636 : : !c data minmag(1),minmag(2) / 128, 0 /
637 : : !c data minmag(3),minmag(4) / 0, 0 /
638 : : !c
639 : : !c data maxmag(1),maxmag(2) / 32767, -1 /
640 : : !c data maxmag(3),maxmag(4) / -1, -1 /
641 : : !c
642 : : !c Machine constants for the Burroughs 6700/7700 systems.
643 : : !c
644 : : !c data mcheps(1) / o1451000000000000 /
645 : : !c data mcheps(2) / o0000000000000000 /
646 : : !c
647 : : !c data minmag(1) / o1771000000000000 /
648 : : !c data minmag(2) / o7770000000000000 /
649 : : !c
650 : : !c data maxmag(1) / o0777777777777777 /
651 : : !c data maxmag(2) / o7777777777777777 /
652 : : !c
653 : : !c Machine constants for the Burroughs 5700 system.
654 : : !c
655 : : !c data mcheps(1) / o1451000000000000 /
656 : : !c data mcheps(2) / o0000000000000000 /
657 : : !c
658 : : !c data minmag(1) / o1771000000000000 /
659 : : !c data minmag(2) / o0000000000000000 /
660 : : !c
661 : : !c data maxmag(1) / o0777777777777777 /
662 : : !c data maxmag(2) / o0007777777777777 /
663 : : !c
664 : : !c Machine constants for the Burroughs 1700 system.
665 : : !c
666 : : !c data mcheps(1) / zcc6800000 /
667 : : !c data mcheps(2) / z000000000 /
668 : : !c
669 : : !c data minmag(1) / zc00800000 /
670 : : !c data minmag(2) / z000000000 /
671 : : !c
672 : : !c data maxmag(1) / zdffffffff /
673 : : !c data maxmag(2) / zfffffffff /
674 : : !c
675 : : !c Machine constants for the Univac 1100 series.
676 : : !c
677 : : !c data mcheps(1),mcheps(2) / o170640000000, o000000000000 /
678 : : !c data minmag(1),minmag(2) / o000040000000, o000000000000 /
679 : : !c data maxmag(1),maxmag(2) / o377777777777, o777777777777 /
680 : : !c
681 : : !c Machine constants for the Data General Eclipse S/200.
682 : : !c
683 : : !c Note - it may be appropriate to include the following card -
684 : : !c static dmach(3)
685 : : !c
686 : : !c data minmag/20k,3*0/,maxmag/77777k,3*177777k/
687 : : !c data mcheps/32020k,3*0/
688 : : !c
689 : : !c Machine constants for the Harris 220.
690 : : !c
691 : : !c data mcheps(1),mcheps(2) / '20000000, '00000334 /
692 : : !c data minmag(1),minmag(2) / '20000000, '00000201 /
693 : : !c data maxmag(1),maxmag(2) / '37777777, '37777577 /
694 : : !c
695 : : !c Machine constants for the Cray-1.
696 : : !c
697 : : !c data mcheps(1) / 0376424000000000000000b /
698 : : !c data mcheps(2) / 0000000000000000000000b /
699 : : !c
700 : : !c data minmag(1) / 0200034000000000000000b /
701 : : !c data minmag(2) / 0000000000000000000000b /
702 : : !c
703 : : !c data maxmag(1) / 0577777777777777777777b /
704 : : !c data maxmag(2) / 0000007777777777777776b /
705 : : !c
706 : : !c Machine constants for the Prime 400.
707 : : !c
708 : : !c data mcheps(1),mcheps(2) / :10000000000, :00000000123 /
709 : : !c data minmag(1),minmag(2) / :10000000000, :00000100000 /
710 : : !c data maxmag(1),maxmag(2) / :17777777777, :37777677776 /
711 : : !c
712 : : !c Machine constants for the VAX-11.
713 : : !c
714 : : !c data mcheps(1),mcheps(2) / 9472, 0 /
715 : : !c data minmag(1),minmag(2) / 128, 0 /
716 : : !c data maxmag(1),maxmag(2) / -32769, -1 /
717 : : !c
718 : : !c Machine constants for IEEE machines.
719 : : !c
720 : : data dmach(1) /2.22044604926d-16/
721 : : data dmach(2) /2.22507385852d-308/
722 : : data dmach(3) /1.79769313485d+308/
723 : : !c
724 [ # # ]: 0 : dpmpar = dmach(i)
725 : : return
726 : : !c
727 : : !c Last card of function dpmpar.
728 : : !c
729 : : end function dpmpar
730 : :
731 : : recursive &
732 : 0 : &double precision function enorm(n,x)
733 : : integer n
734 : : double precision x(n)
735 : : !c **********
736 : : !c
737 : : !c function enorm
738 : : !c
739 : : !c given an n-vector x, this function calculates the
740 : : !c euclidean norm of x.
741 : : !c
742 : : !c the euclidean norm is computed by accumulating the sum of
743 : : !c squares in three different sums. the sums of squares for the
744 : : !c small and large components are scaled so that no overflows
745 : : !c occur. non-destructive underflows are permitted. underflows
746 : : !c and overflows do not occur in the computation of the unscaled
747 : : !c sum of squares for the intermediate components.
748 : : !c the definitions of small, intermediate and large components
749 : : !c depend on two constants, rdwarf and rgiant. the main
750 : : !c restrictions on these constants are that rdwarf**2 not
751 : : !c underflow and rgiant**2 not overflow. the constants
752 : : !c given here are suitable for every known computer.
753 : : !c
754 : : !c the function statement is
755 : : !c
756 : : !c double precision function enorm(n,x)
757 : : !c
758 : : !c where
759 : : !c
760 : : !c n is a positive integer input variable.
761 : : !c
762 : : !c x is an input array of length n.
763 : : !c
764 : : !c subprograms called
765 : : !c
766 : : !c fortran-supplied ... dabs,dsqrt
767 : : !c
768 : : !c argonne national laboratory. minpack project. march 1980.
769 : : !c burton s. garbow, kenneth e. hillstrom, jorge j. more
770 : : !c
771 : : !c **********
772 : : integer i
773 : : double precision agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs, &
774 : : & x1max,x3max,zero
775 : : data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/
776 : 0 : s1 = zero
777 : 0 : s2 = zero
778 : 0 : s3 = zero
779 : 0 : x1max = zero
780 : 0 : x3max = zero
781 : 0 : floatn = n
782 : 0 : agiant = rgiant/floatn
783 [ # # ]: 0 : do 90 i = 1, n
784 : 0 : xabs = dabs(x(i))
785 [ # # # # ]: 0 : if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70
786 [ # # ]: 0 : if (xabs .le. rdwarf) go to 30
787 : : !c
788 : : !c sum for large components.
789 : : !c
790 [ # # ]: 0 : if (xabs .le. x1max) go to 10
791 : 0 : s1 = one + s1*(x1max/xabs)**2
792 : 0 : x1max = xabs
793 : 0 : go to 20
794 : : 10 continue
795 : 0 : s1 = s1 + (xabs/x1max)**2
796 : : 20 continue
797 : 0 : go to 60
798 : : 30 continue
799 : : !c
800 : : !c sum for small components.
801 : : !c
802 [ # # ]: 0 : if (xabs .le. x3max) go to 40
803 : 0 : s3 = one + s3*(x3max/xabs)**2
804 : 0 : x3max = xabs
805 : 0 : go to 50
806 : : 40 continue
807 [ # # ]: 0 : if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2
808 : : 50 continue
809 : : 60 continue
810 : 0 : go to 80
811 : : 70 continue
812 : : !c
813 : : !c sum for intermediate components.
814 : : !c
815 : 0 : s2 = s2 + xabs**2
816 : : 80 continue
817 : 0 : 90 continue
818 : : !c
819 : : !c calculation of norm.
820 : : !c
821 [ # # ]: 0 : if (s1 .eq. zero) go to 100
822 : 0 : enorm = x1max*dsqrt(s1+(s2/x1max)/x1max)
823 : 0 : go to 130
824 : : 100 continue
825 [ # # ]: 0 : if (s2 .eq. zero) go to 110
826 [ # # ]: 0 : if (s2 .ge. x3max) &
827 : 0 : & enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3)))
828 [ # # ]: 0 : if (s2 .lt. x3max) &
829 : 0 : & enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3)))
830 : 0 : go to 120
831 : : 110 continue
832 : 0 : enorm = x3max*dsqrt(s3)
833 : : 120 continue
834 : : 130 continue
835 : : return
836 : : !c
837 : : !c last card of function enorm.
838 : : !c
839 : : end function enorm
840 : :
841 : : recursive &
842 : 0 : &subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn, &
843 : 0 : & wa1,wa2)
844 : : integer n,ldfjac,iflag,ml,mu
845 : : double precision epsfcn
846 : : double precision x(n),fvec(n),fjac(ldfjac,n),wa1(n),wa2(n)
847 : : interface
848 : : subroutine fcn(n,x,fvec,iflag)
849 : : integer, intent(in) :: n,iflag
850 : : double precision, intent(in) :: x(n)
851 : : double precision, intent(inout) :: fvec(n)
852 : : end subroutine fcn
853 : : end interface
854 : : !c **********
855 : : !c
856 : : !c subroutine fdjac1
857 : : !c
858 : : !c this subroutine computes a forward-difference approximation
859 : : !c to the n by n jacobian matrix associated with a specified
860 : : !c problem of n functions in n variables. if the jacobian has
861 : : !c a banded form, then function evaluations are saved by only
862 : : !c approximating the nonzero terms.
863 : : !c
864 : : !c the subroutine statement is
865 : : !c
866 : : !c subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
867 : : !c wa1,wa2)
868 : : !c
869 : : !c where
870 : : !c
871 : : !c fcn is the name of the user-supplied subroutine which
872 : : !c calculates the functions. fcn must be declared
873 : : !c in an external statement in the user calling
874 : : !c program, and should be written as follows.
875 : : !c
876 : : !c subroutine fcn(n,x,fvec,iflag)
877 : : !c integer n,iflag
878 : : !c double precision x(n),fvec(n)
879 : : !c ----------
880 : : !c calculate the functions at x and
881 : : !c return this vector in fvec.
882 : : !c ----------
883 : : !c return
884 : : !c end
885 : : !c
886 : : !c the value of iflag should not be changed by fcn unless
887 : : !c the user wants to terminate execution of fdjac1.
888 : : !c in this case set iflag to a negative integer.
889 : : !c
890 : : !c n is a positive integer input variable set to the number
891 : : !c of functions and variables.
892 : : !c
893 : : !c x is an input array of length n.
894 : : !c
895 : : !c fvec is an input array of length n which must contain the
896 : : !c functions evaluated at x.
897 : : !c
898 : : !c fjac is an output n by n array which contains the
899 : : !c approximation to the jacobian matrix evaluated at x.
900 : : !c
901 : : !c ldfjac is a positive integer input variable not less than n
902 : : !c which specifies the leading dimension of the array fjac.
903 : : !c
904 : : !c iflag is an integer variable which can be used to terminate
905 : : !c the execution of fdjac1. see description of fcn.
906 : : !c
907 : : !c ml is a nonnegative integer input variable which specifies
908 : : !c the number of subdiagonals within the band of the
909 : : !c jacobian matrix. if the jacobian is not banded, set
910 : : !c ml to at least n - 1.
911 : : !c
912 : : !c epsfcn is an input variable used in determining a suitable
913 : : !c step length for the forward-difference approximation. this
914 : : !c approximation assumes that the relative errors in the
915 : : !c functions are of the order of epsfcn. if epsfcn is less
916 : : !c than the machine precision, it is assumed that the relative
917 : : !c errors in the functions are of the order of the machine
918 : : !c precision.
919 : : !c
920 : : !c mu is a nonnegative integer input variable which specifies
921 : : !c the number of superdiagonals within the band of the
922 : : !c jacobian matrix. if the jacobian is not banded, set
923 : : !c mu to at least n - 1.
924 : : !c
925 : : !c wa1 and wa2 are work arrays of length n. if ml + mu + 1 is at
926 : : !c least n, then the jacobian is considered dense, and wa2 is
927 : : !c not referenced.
928 : : !c
929 : : !c subprograms called
930 : : !c
931 : : !c minpack-supplied ... dpmpar
932 : : !c
933 : : !c fortran-supplied ... dabs,dmax1,dsqrt
934 : : !c
935 : : !c argonne national laboratory. minpack project. march 1980.
936 : : !c burton s. garbow, kenneth e. hillstrom, jorge j. more
937 : : !c
938 : : !c **********
939 : : integer i,j,k,msum
940 : : double precision eps,epsmch,h,temp,zero
941 : : !double precision dpmpar
942 : : data zero /0.0d0/
943 : : !c
944 : : !c epsmch is the machine precision.
945 : : !c
946 : 0 : epsmch = dpmpar(1)
947 : : !c
948 : 0 : eps = dsqrt(dmax1(epsfcn,epsmch))
949 : 0 : msum = ml + mu + 1
950 [ # # ]: 0 : if (msum .lt. n) go to 40
951 : : !c
952 : : !c computation of dense approximate jacobian.
953 : : !c
954 [ # # ]: 0 : do 20 j = 1, n
955 : 0 : temp = x(j)
956 : 0 : h = eps*dabs(temp)
957 [ # # ]: 0 : if (h .eq. zero) h = eps
958 : 0 : x(j) = temp + h
959 : 0 : call fcn(n,x,wa1,iflag)
960 [ # # ]: 0 : if (iflag .lt. 0) go to 30
961 : 0 : x(j) = temp
962 [ # # ]: 0 : do 10 i = 1, n
963 : 0 : fjac(i,j) = (wa1(i) - fvec(i))/h
964 : 0 : 10 continue
965 : 0 : 20 continue
966 : : 30 continue
967 : 0 : go to 110
968 : : 40 continue
969 : : !c
970 : : !c computation of banded approximate jacobian.
971 : : !c
972 [ # # ]: 0 : do 90 k = 1, msum
973 [ # # ]: 0 : do 60 j = k, n, msum
974 : 0 : wa2(j) = x(j)
975 : 0 : h = eps*dabs(wa2(j))
976 [ # # ]: 0 : if (h .eq. zero) h = eps
977 [ # # ]: 0 : x(j) = wa2(j) + h
978 : 0 : 60 continue
979 : 0 : call fcn(n,x,wa1,iflag)
980 [ # # ]: 0 : if (iflag .lt. 0) go to 100
981 [ # # # # ]: 0 : do 80 j = k, n, msum
982 : 0 : x(j) = wa2(j)
983 : 0 : h = eps*dabs(wa2(j))
984 [ # # ]: 0 : if (h .eq. zero) h = eps
985 [ # # ]: 0 : do 70 i = 1, n
986 : 0 : fjac(i,j) = zero
987 [ # # # # ]: 0 : if (i .ge. j - mu .and. i .le. j + ml) &
988 : 0 : & fjac(i,j) = (wa1(i) - fvec(i))/h
989 : 0 : 70 continue
990 : 0 : 80 continue
991 : 0 : 90 continue
992 : : 100 continue
993 : : 110 continue
994 : 0 : return
995 : : !c
996 : : !c last card of subroutine fdjac1.
997 : : !c
998 : : end subroutine fdjac1
999 : :
1000 : : recursive &
1001 : 0 : &subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
1002 : : integer m,n,lda,lipvt
1003 : : integer ipvt(lipvt)
1004 : : logical pivot
1005 : : double precision a(lda,n),rdiag(n),acnorm(n),wa(n)
1006 : : !c **********
1007 : : !c
1008 : : !c subroutine qrfac
1009 : : !c
1010 : : !c this subroutine uses householder transformations with column
1011 : : !c pivoting (optional) to compute a qr factorization of the
1012 : : !c m by n matrix a. that is, qrfac determines an orthogonal
1013 : : !c matrix q, a permutation matrix p, and an upper trapezoidal
1014 : : !c matrix r with diagonal elements of nonincreasing magnitude,
1015 : : !c such that a*p = q*r. the householder transformation for
1016 : : !c column k, k = 1,2,...,min(m,n), is of the form
1017 : : !c
1018 : : !c t
1019 : : !c i - (1/u(k))*u*u
1020 : : !c
1021 : : !c where u has zeros in the first k-1 positions. the form of
1022 : : !c this transformation and the method of pivoting first
1023 : : !c appeared in the corresponding linpack subroutine.
1024 : : !c
1025 : : !c the subroutine statement is
1026 : : !c
1027 : : !c subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
1028 : : !c
1029 : : !c where
1030 : : !c
1031 : : !c m is a positive integer input variable set to the number
1032 : : !c of rows of a.
1033 : : !c
1034 : : !c n is a positive integer input variable set to the number
1035 : : !c of columns of a.
1036 : : !c
1037 : : !c a is an m by n array. on input a contains the matrix for
1038 : : !c which the qr factorization is to be computed. on output
1039 : : !c the strict upper trapezoidal part of a contains the strict
1040 : : !c upper trapezoidal part of r, and the lower trapezoidal
1041 : : !c part of a contains a factored form of q (the non-trivial
1042 : : !c elements of the u vectors described above).
1043 : : !c
1044 : : !c lda is a positive integer input variable not less than m
1045 : : !c which specifies the leading dimension of the array a.
1046 : : !c
1047 : : !c pivot is a logical input variable. if pivot is set true,
1048 : : !c then column pivoting is enforced. if pivot is set false,
1049 : : !c then no column pivoting is done.
1050 : : !c
1051 : : !c ipvt is an integer output array of length lipvt. ipvt
1052 : : !c defines the permutation matrix p such that a*p = q*r.
1053 : : !c column j of p is column ipvt(j) of the identity matrix.
1054 : : !c if pivot is false, ipvt is not referenced.
1055 : : !c
1056 : : !c lipvt is a positive integer input variable. if pivot is false,
1057 : : !c then lipvt may be as small as 1. if pivot is true, then
1058 : : !c lipvt must be at least n.
1059 : : !c
1060 : : !c rdiag is an output array of length n which contains the
1061 : : !c diagonal elements of r.
1062 : : !c
1063 : : !c acnorm is an output array of length n which contains the
1064 : : !c norms of the corresponding columns of the input matrix a.
1065 : : !c if this information is not needed, then acnorm can coincide
1066 : : !c with rdiag.
1067 : : !c
1068 : : !c wa is a work array of length n. if pivot is false, then wa
1069 : : !c can coincide with rdiag.
1070 : : !c
1071 : : !c subprograms called
1072 : : !c
1073 : : !c minpack-supplied ... dpmpar,enorm
1074 : : !c
1075 : : !c fortran-supplied ... dmax1,dsqrt,min0
1076 : : !c
1077 : : !c argonne national laboratory. minpack project. march 1980.
1078 : : !c burton s. garbow, kenneth e. hillstrom, jorge j. more
1079 : : !c
1080 : : !c **********
1081 : : integer i,j,jp1,k,kmax,minmn
1082 : : double precision ajnorm,epsmch,one,p05,sum,temp,zero
1083 : : !double precision dpmpar,enorm
1084 : : data one,p05,zero /1.0d0,5.0d-2,0.0d0/
1085 : : !c
1086 : : !c epsmch is the machine precision.
1087 : : !c
1088 : 0 : epsmch = dpmpar(1)
1089 : : !c
1090 : : !c compute the initial column norms and initialize several arrays.
1091 : : !c
1092 [ # # ]: 0 : do 10 j = 1, n
1093 : 0 : acnorm(j) = enorm(m,a(1,j))
1094 : 0 : rdiag(j) = acnorm(j)
1095 : 0 : wa(j) = rdiag(j)
1096 [ # # ]: 0 : if (pivot) ipvt(j) = j
1097 : 0 : 10 continue
1098 : : !c
1099 : : !c reduce a to r with householder transformations.
1100 : : !c
1101 : 0 : minmn = min0(m,n)
1102 [ # # ]: 0 : do 110 j = 1, minmn
1103 [ # # ]: 0 : if (.not.pivot) go to 40
1104 : : !c
1105 : : !c bring the column of largest norm into the pivot position.
1106 : : !c
1107 : : kmax = j
1108 [ # # ]: 0 : do 20 k = j, n
1109 [ # # ]: 0 : if (rdiag(k) .gt. rdiag(kmax)) kmax = k
1110 : 0 : 20 continue
1111 [ # # ]: 0 : if (kmax .eq. j) go to 40
1112 [ # # ]: 0 : do 30 i = 1, m
1113 : 0 : temp = a(i,j)
1114 : 0 : a(i,j) = a(i,kmax)
1115 : 0 : a(i,kmax) = temp
1116 : 0 : 30 continue
1117 : 0 : rdiag(kmax) = rdiag(j)
1118 : 0 : wa(kmax) = wa(j)
1119 : 0 : k = ipvt(j)
1120 : 0 : ipvt(j) = ipvt(kmax)
1121 : 0 : ipvt(kmax) = k
1122 : : 40 continue
1123 : : !c
1124 : : !c compute the householder transformation to reduce the
1125 : : !c j-th column of a to a multiple of the j-th unit vector.
1126 : : !c
1127 : 0 : ajnorm = enorm(m-j+1,a(j,j))
1128 [ # # ]: 0 : if (ajnorm .eq. zero) go to 100
1129 [ # # ]: 0 : if (a(j,j) .lt. zero) ajnorm = -ajnorm
1130 [ # # ]: 0 : do 50 i = j, m
1131 : 0 : a(i,j) = a(i,j)/ajnorm
1132 : 0 : 50 continue
1133 : 0 : a(j,j) = a(j,j) + one
1134 : : !c
1135 : : !c apply the transformation to the remaining columns
1136 : : !c and update the norms.
1137 : : !c
1138 : 0 : jp1 = j + 1
1139 [ # # ]: 0 : if (n .lt. jp1) go to 100
1140 [ # # ]: 0 : do 90 k = jp1, n
1141 : : sum = zero
1142 [ # # ]: 0 : do 60 i = j, m
1143 : 0 : sum = sum + a(i,j)*a(i,k)
1144 : 0 : 60 continue
1145 : 0 : temp = sum/a(j,j)
1146 [ # # ]: 0 : do 70 i = j, m
1147 : 0 : a(i,k) = a(i,k) - temp*a(i,j)
1148 : 0 : 70 continue
1149 [ # # # # ]: 0 : if (.not.pivot .or. rdiag(k) .eq. zero) go to 80
1150 : 0 : temp = a(j,k)/rdiag(k)
1151 : 0 : rdiag(k) = rdiag(k)*dsqrt(dmax1(zero,one-temp**2))
1152 [ # # ]: 0 : if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80
1153 : 0 : rdiag(k) = enorm(m-j,a(jp1,k))
1154 : 0 : wa(k) = rdiag(k)
1155 : : 80 continue
1156 : 0 : 90 continue
1157 : : 100 continue
1158 : 0 : rdiag(j) = -ajnorm
1159 : 0 : 110 continue
1160 : 0 : return
1161 : : !c
1162 : : !c last card of subroutine qrfac.
1163 : : !c
1164 : : end subroutine qrfac
1165 : :
1166 : : recursive &
1167 : 0 : &subroutine qform(m,n,q,ldq,wa)
1168 : : integer m,n,ldq
1169 : : double precision q(ldq,m),wa(m)
1170 : : !c **********
1171 : : !c
1172 : : !c subroutine qform
1173 : : !c
1174 : : !c this subroutine proceeds from the computed qr factorization of
1175 : : !c an m by n matrix a to accumulate the m by m orthogonal matrix
1176 : : !c q from its factored form.
1177 : : !c
1178 : : !c the subroutine statement is
1179 : : !c
1180 : : !c subroutine qform(m,n,q,ldq,wa)
1181 : : !c
1182 : : !c where
1183 : : !c
1184 : : !c m is a positive integer input variable set to the number
1185 : : !c of rows of a and the order of q.
1186 : : !c
1187 : : !c n is a positive integer input variable set to the number
1188 : : !c of columns of a.
1189 : : !c
1190 : : !c q is an m by m array. on input the full lower trapezoid in
1191 : : !c the first min(m,n) columns of q contains the factored form.
1192 : : !c on output q has been accumulated into a square matrix.
1193 : : !c
1194 : : !c ldq is a positive integer input variable not less than m
1195 : : !c which specifies the leading dimension of the array q.
1196 : : !c
1197 : : !c wa is a work array of length m.
1198 : : !c
1199 : : !c subprograms called
1200 : : !c
1201 : : !c fortran-supplied ... min0
1202 : : !c
1203 : : !c argonne national laboratory. minpack project. march 1980.
1204 : : !c burton s. garbow, kenneth e. hillstrom, jorge j. more
1205 : : !c
1206 : : !c **********
1207 : : integer i,j,jm1,k,l,minmn,np1
1208 : : double precision one,sum,temp,zero
1209 : : data one,zero /1.0d0,0.0d0/
1210 : : !c
1211 : : !c zero out upper triangle of q in the first min(m,n) columns.
1212 : : !c
1213 : 0 : minmn = min0(m,n)
1214 [ # # ]: 0 : if (minmn .lt. 2) go to 30
1215 [ # # ]: 0 : do 20 j = 2, minmn
1216 : : jm1 = j - 1
1217 [ # # ]: 0 : do 10 i = 1, jm1
1218 : 0 : q(i,j) = zero
1219 : 0 : 10 continue
1220 : 0 : 20 continue
1221 : : 30 continue
1222 : : !c
1223 : : !c initialize remaining columns to those of the identity matrix.
1224 : : !c
1225 : 0 : np1 = n + 1
1226 [ # # ]: 0 : if (m .lt. np1) go to 60
1227 [ # # ]: 0 : do 50 j = np1, m
1228 [ # # ]: 0 : do 40 i = 1, m
1229 : 0 : q(i,j) = zero
1230 : 0 : 40 continue
1231 : 0 : q(j,j) = one
1232 : 0 : 50 continue
1233 : : 60 continue
1234 : : !c
1235 : : !c accumulate q from its factored form.
1236 : : !c
1237 [ # # ]: 0 : do 120 l = 1, minmn
1238 : 0 : k = minmn - l + 1
1239 [ # # ]: 0 : do 70 i = k, m
1240 : 0 : wa(i) = q(i,k)
1241 : 0 : q(i,k) = zero
1242 : 0 : 70 continue
1243 : 0 : q(k,k) = one
1244 [ # # ]: 0 : if (wa(k) .eq. zero) go to 110
1245 [ # # ]: 0 : do 100 j = k, m
1246 : : sum = zero
1247 [ # # ]: 0 : do 80 i = k, m
1248 : 0 : sum = sum + q(i,j)*wa(i)
1249 : 0 : 80 continue
1250 : 0 : temp = sum/wa(k)
1251 [ # # ]: 0 : do 90 i = k, m
1252 : 0 : q(i,j) = q(i,j) - temp*wa(i)
1253 : 0 : 90 continue
1254 : 0 : 100 continue
1255 : : 110 continue
1256 : 0 : 120 continue
1257 : 0 : return
1258 : : !c
1259 : : !c last card of subroutine qform.
1260 : : !c
1261 : : end subroutine qform
1262 : :
1263 : : recursive &
1264 : 0 : &subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
1265 : : integer n,lr
1266 : : double precision delta
1267 : : double precision r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n)
1268 : : !c **********
1269 : : !c
1270 : : !c subroutine dogleg
1271 : : !c
1272 : : !c given an m by n matrix a, an n by n nonsingular diagonal
1273 : : !c matrix d, an m-vector b, and a positive number delta, the
1274 : : !c problem is to determine the convex combination x of the
1275 : : !c gauss-newton and scaled gradient directions that minimizes
1276 : : !c (a*x - b) in the least squares sense, subject to the
1277 : : !c restriction that the euclidean norm of d*x be at most delta.
1278 : : !c
1279 : : !c this subroutine completes the solution of the problem
1280 : : !c if it is provided with the necessary information from the
1281 : : !c qr factorization of a. that is, if a = q*r, where q has
1282 : : !c orthogonal columns and r is an upper triangular matrix,
1283 : : !c then dogleg expects the full upper triangle of r and
1284 : : !c the first n components of (q transpose)*b.
1285 : : !c
1286 : : !c the subroutine statement is
1287 : : !c
1288 : : !c subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
1289 : : !c
1290 : : !c where
1291 : : !c
1292 : : !c n is a positive integer input variable set to the order of r.
1293 : : !c
1294 : : !c r is an input array of length lr which must contain the upper
1295 : : !c triangular matrix r stored by rows.
1296 : : !c
1297 : : !c lr is a positive integer input variable not less than
1298 : : !c (n*(n+1))/2.
1299 : : !c
1300 : : !c diag is an input array of length n which must contain the
1301 : : !c diagonal elements of the matrix d.
1302 : : !c
1303 : : !c qtb is an input array of length n which must contain the first
1304 : : !c n elements of the vector (q transpose)*b.
1305 : : !c
1306 : : !c delta is a positive input variable which specifies an upper
1307 : : !c bound on the euclidean norm of d*x.
1308 : : !c
1309 : : !c x is an output array of length n which contains the desired
1310 : : !c convex combination of the gauss-newton direction and the
1311 : : !c scaled gradient direction.
1312 : : !c
1313 : : !c wa1 and wa2 are work arrays of length n.
1314 : : !c
1315 : : !c subprograms called
1316 : : !c
1317 : : !c minpack-supplied ... dpmpar,enorm
1318 : : !c
1319 : : !c fortran-supplied ... dabs,dmax1,dmin1,dsqrt
1320 : : !c
1321 : : !c argonne national laboratory. minpack project. march 1980.
1322 : : !c burton s. garbow, kenneth e. hillstrom, jorge j. more
1323 : : !c
1324 : : !c **********
1325 : : integer i,j,jj,jp1,k,l
1326 : : double precision alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum, &
1327 : : & temp,zero
1328 : : !double precision dpmpar,enorm
1329 : : data one,zero /1.0d0,0.0d0/
1330 : : !c
1331 : : !c epsmch is the machine precision.
1332 : : !c
1333 : 0 : epsmch = dpmpar(1)
1334 : : !c
1335 : : !c first, calculate the gauss-newton direction.
1336 : : !c
1337 : 0 : jj = (n*(n + 1))/2 + 1
1338 [ # # ]: 0 : do 50 k = 1, n
1339 : 0 : j = n - k + 1
1340 : 0 : jp1 = j + 1
1341 : 0 : jj = jj - k
1342 : 0 : l = jj + 1
1343 : 0 : sum = zero
1344 [ # # ]: 0 : if (n .lt. jp1) go to 20
1345 [ # # ]: 0 : do 10 i = jp1, n
1346 : 0 : sum = sum + r(l)*x(i)
1347 : 0 : l = l + 1
1348 : 0 : 10 continue
1349 : : 20 continue
1350 : 0 : temp = r(jj)
1351 [ # # ]: 0 : if (temp .ne. zero) go to 40
1352 : : l = j
1353 [ # # ]: 0 : do 30 i = 1, j
1354 : 0 : temp = dmax1(temp,dabs(r(l)))
1355 : 0 : l = l + n - i
1356 : 0 : 30 continue
1357 : 0 : temp = epsmch*temp
1358 [ # # ]: 0 : if (temp .eq. zero) temp = epsmch
1359 : : 40 continue
1360 : 0 : x(j) = (qtb(j) - sum)/temp
1361 : 0 : 50 continue
1362 : : !c
1363 : : !c test whether the gauss-newton direction is acceptable.
1364 : : !c
1365 [ # # ]: 0 : do 60 j = 1, n
1366 : 0 : wa1(j) = zero
1367 : 0 : wa2(j) = diag(j)*x(j)
1368 : 0 : 60 continue
1369 : 0 : qnorm = enorm(n,wa2)
1370 [ # # ]: 0 : if (qnorm .le. delta) go to 140
1371 : : !c
1372 : : !c the gauss-newton direction is not acceptable.
1373 : : !c next, calculate the scaled gradient direction.
1374 : : !c
1375 : : l = 1
1376 [ # # ]: 0 : do 80 j = 1, n
1377 : 0 : temp = qtb(j)
1378 [ # # ]: 0 : do 70 i = j, n
1379 : 0 : wa1(i) = wa1(i) + r(l)*temp
1380 : 0 : l = l + 1
1381 : 0 : 70 continue
1382 : 0 : wa1(j) = wa1(j)/diag(j)
1383 : 0 : 80 continue
1384 : : !c
1385 : : !c calculate the norm of the scaled gradient and test for
1386 : : !c the special case in which the scaled gradient is zero.
1387 : : !c
1388 : 0 : gnorm = enorm(n,wa1)
1389 : 0 : sgnorm = zero
1390 : 0 : alpha = delta/qnorm
1391 [ # # ]: 0 : if (gnorm .eq. zero) go to 120
1392 : : !c
1393 : : !c calculate the point along the scaled gradient
1394 : : !c at which the quadratic is minimized.
1395 : : !c
1396 [ # # ]: 0 : do 90 j = 1, n
1397 : 0 : wa1(j) = (wa1(j)/gnorm)/diag(j)
1398 : 0 : 90 continue
1399 : : l = 1
1400 [ # # ]: 0 : do 110 j = 1, n
1401 : : sum = zero
1402 [ # # ]: 0 : do 100 i = j, n
1403 : 0 : sum = sum + r(l)*wa1(i)
1404 : 0 : l = l + 1
1405 : 0 : 100 continue
1406 : 0 : wa2(j) = sum
1407 : 0 : 110 continue
1408 : 0 : temp = enorm(n,wa2)
1409 : 0 : sgnorm = (gnorm/temp)/temp
1410 : : !c
1411 : : !c test whether the scaled gradient direction is acceptable.
1412 : : !c
1413 : 0 : alpha = zero
1414 [ # # ]: 0 : if (sgnorm .ge. delta) go to 120
1415 : : !c
1416 : : !c the scaled gradient direction is not acceptable.
1417 : : !c finally, calculate the point along the dogleg
1418 : : !c at which the quadratic is minimized.
1419 : : !c
1420 : 0 : bnorm = enorm(n,qtb)
1421 : 0 : temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
1422 : : temp = temp - (delta/qnorm)*(sgnorm/delta)**2 &
1423 : : & + dsqrt((temp-(delta/qnorm))**2 &
1424 : 0 : & +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2))
1425 : 0 : alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp
1426 : : 120 continue
1427 : : !c
1428 : : !c form appropriate convex combination of the gauss-newton
1429 : : !c direction and the scaled gradient direction.
1430 : : !c
1431 : 0 : temp = (one - alpha)*dmin1(sgnorm,delta)
1432 [ # # ]: 0 : do 130 j = 1, n
1433 : 0 : x(j) = temp*wa1(j) + alpha*x(j)
1434 : 0 : 130 continue
1435 : : 140 continue
1436 : 0 : return
1437 : : !c
1438 : : !c last card of subroutine dogleg.
1439 : : !c
1440 : : end subroutine dogleg
1441 : :
1442 : : recursive &
1443 : 0 : &subroutine r1updt(m,n,s,ls,u,v,w,sing)
1444 : : integer m,n,ls
1445 : : logical sing
1446 : : double precision s(ls),u(m),v(n),w(m)
1447 : : !c **********
1448 : : !c
1449 : : !c subroutine r1updt
1450 : : !c
1451 : : !c given an m by n lower trapezoidal matrix s, an m-vector u,
1452 : : !c and an n-vector v, the problem is to determine an
1453 : : !c orthogonal matrix q such that
1454 : : !c
1455 : : !c t
1456 : : !c (s + u*v )*q
1457 : : !c
1458 : : !c is again lower trapezoidal.
1459 : : !c
1460 : : !c this subroutine determines q as the product of 2*(n - 1)
1461 : : !c transformations
1462 : : !c
1463 : : !c gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
1464 : : !c
1465 : : !c where gv(i), gw(i) are givens rotations in the (i,n) plane
1466 : : !c which eliminate elements in the i-th and n-th planes,
1467 : : !c respectively. q itself is not accumulated, rather the
1468 : : !c information to recover the gv, gw rotations is returned.
1469 : : !c
1470 : : !c the subroutine statement is
1471 : : !c
1472 : : !c subroutine r1updt(m,n,s,ls,u,v,w,sing)
1473 : : !c
1474 : : !c where
1475 : : !c
1476 : : !c m is a positive integer input variable set to the number
1477 : : !c of rows of s.
1478 : : !c
1479 : : !c n is a positive integer input variable set to the number
1480 : : !c of columns of s. n must not exceed m.
1481 : : !c
1482 : : !c s is an array of length ls. on input s must contain the lower
1483 : : !c trapezoidal matrix s stored by columns. on output s contains
1484 : : !c the lower trapezoidal matrix produced as described above.
1485 : : !c
1486 : : !c ls is a positive integer input variable not less than
1487 : : !c (n*(2*m-n+1))/2.
1488 : : !c
1489 : : !c u is an input array of length m which must contain the
1490 : : !c vector u.
1491 : : !c
1492 : : !c v is an array of length n. on input v must contain the vector!
1493 : : !c v. on output v(i) contains the information necessary to
1494 : : !c recover the givens rotation gv(i) described above.
1495 : : !c
1496 : : !c w is an output array of length m. w(i) contains information
1497 : : !c necessary to recover the givens rotation gw(i) described
1498 : : !c above.
1499 : : !c
1500 : : !c sing is a logical output variable. sing is set true if any
1501 : : !c of the diagonal elements of the output s are zero. otherwise
1502 : : !c sing is set false.
1503 : : !c
1504 : : !c subprograms called
1505 : : !c
1506 : : !c minpack-supplied ... dpmpar
1507 : : !c
1508 : : !c fortran-supplied ... dabs,dsqrt
1509 : : !c
1510 : : !c argonne national laboratory. minpack project. march 1980.
1511 : : !c burton s. garbow, kenneth e. hillstrom, jorge j. more,
1512 : : !c john l. nazareth
1513 : : !c
1514 : : !c **********
1515 : : integer i,j,jj,l,nmj,nm1
1516 : : double precision cos,cotan,giant,one,p5,p25,sin,tan,tau,temp, &
1517 : : & zero
1518 : : !double precision dpmpar
1519 : : data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
1520 : : !c
1521 : : !c giant is the largest magnitude.
1522 : : !c
1523 : 0 : giant = dpmpar(3)
1524 : : !c
1525 : : !c initialize the diagonal element pointer.
1526 : : !c
1527 : 0 : jj = (n*(2*m - n + 1))/2 - (m - n)
1528 : : !c
1529 : : !c move the nontrivial part of the last column of s into w.
1530 : : !c
1531 : 0 : l = jj
1532 [ # # ]: 0 : do 10 i = n, m
1533 : 0 : w(i) = s(l)
1534 : 0 : l = l + 1
1535 : 0 : 10 continue
1536 : : !c
1537 : : !c rotate the vector v into a multiple of the n-th unit vector
1538 : : !c in such a way that a spike is introduced into w.
1539 : : !c
1540 : 0 : nm1 = n - 1
1541 [ # # ]: 0 : if (nm1 .lt. 1) go to 70
1542 [ # # ]: 0 : do 60 nmj = 1, nm1
1543 : 0 : j = n - nmj
1544 : 0 : jj = jj - (m - j + 1)
1545 : 0 : w(j) = zero
1546 [ # # ]: 0 : if (v(j) .eq. zero) go to 50
1547 : : !c
1548 : : !c determine a givens rotation which eliminates the
1549 : : !c j-th element of v.
1550 : : !c
1551 [ # # ]: 0 : if (dabs(v(n)) .ge. dabs(v(j))) go to 20
1552 : 0 : cotan = v(n)/v(j)
1553 : 0 : sin = p5/dsqrt(p25+p25*cotan**2)
1554 : 0 : cos = sin*cotan
1555 : 0 : tau = one
1556 [ # # ]: 0 : if (dabs(cos)*giant .gt. one) tau = one/cos
1557 : 0 : go to 30
1558 : : 20 continue
1559 : 0 : tan = v(j)/v(n)
1560 : 0 : cos = p5/dsqrt(p25+p25*tan**2)
1561 : 0 : sin = cos*tan
1562 : 0 : tau = sin
1563 : : 30 continue
1564 : : !c
1565 : : !c apply the transformation to v and store the information
1566 : : !c necessary to recover the givens rotation.
1567 : : !c
1568 : 0 : v(n) = sin*v(j) + cos*v(n)
1569 : 0 : v(j) = tau
1570 : : !c
1571 : : !c apply the transformation to s and extend the spike in w.
1572 : : !c
1573 : 0 : l = jj
1574 [ # # ]: 0 : do 40 i = j, m
1575 : 0 : temp = cos*s(l) - sin*w(i)
1576 : 0 : w(i) = sin*s(l) + cos*w(i)
1577 : 0 : s(l) = temp
1578 : 0 : l = l + 1
1579 : 0 : 40 continue
1580 : : 50 continue
1581 : 0 : 60 continue
1582 : : 70 continue
1583 : : !c
1584 : : !c add the spike from the rank 1 update to w.
1585 : : !c
1586 [ # # ]: 0 : do 80 i = 1, m
1587 : 0 : w(i) = w(i) + v(n)*u(i)
1588 : 0 : 80 continue
1589 : : !c
1590 : : !c eliminate the spike.
1591 : : !c
1592 : 0 : sing = .false.
1593 [ # # ]: 0 : if (nm1 .lt. 1) go to 140
1594 [ # # ]: 0 : do 130 j = 1, nm1
1595 [ # # ]: 0 : if (w(j) .eq. zero) go to 120
1596 : : !c
1597 : : !c determine a givens rotation which eliminates the
1598 : : !c j-th element of the spike.
1599 : : !c
1600 [ # # ]: 0 : if (dabs(s(jj)) .ge. dabs(w(j))) go to 90
1601 : 0 : cotan = s(jj)/w(j)
1602 : 0 : sin = p5/dsqrt(p25+p25*cotan**2)
1603 : 0 : cos = sin*cotan
1604 : 0 : tau = one
1605 [ # # ]: 0 : if (dabs(cos)*giant .gt. one) tau = one/cos
1606 : 0 : go to 100
1607 : : 90 continue
1608 : 0 : tan = w(j)/s(jj)
1609 : 0 : cos = p5/dsqrt(p25+p25*tan**2)
1610 : 0 : sin = cos*tan
1611 : 0 : tau = sin
1612 : : 100 continue
1613 : : !c
1614 : : !c apply the transformation to s and reduce the spike in w.
1615 : : !c
1616 : 0 : l = jj
1617 [ # # ]: 0 : do 110 i = j, m
1618 : 0 : temp = cos*s(l) + sin*w(i)
1619 : 0 : w(i) = -sin*s(l) + cos*w(i)
1620 : 0 : s(l) = temp
1621 : 0 : l = l + 1
1622 : 0 : 110 continue
1623 : : !c
1624 : : !c store the information necessary to recover the
1625 : : !c givens rotation.
1626 : : !c
1627 : 0 : w(j) = tau
1628 : : 120 continue
1629 : : !c
1630 : : !c test for zero diagonal elements in the output s.
1631 : : !c
1632 [ # # ]: 0 : if (s(jj) .eq. zero) sing = .true.
1633 : 0 : jj = jj + (m - j + 1)
1634 : 0 : 130 continue
1635 : : 140 continue
1636 : : !c
1637 : : !c move w back into the last column of the output s.
1638 : : !c
1639 : 0 : l = jj
1640 [ # # ]: 0 : do 150 i = n, m
1641 : 0 : s(l) = w(i)
1642 : 0 : l = l + 1
1643 : 0 : 150 continue
1644 [ # # ]: 0 : if (s(jj) .eq. zero) sing = .true.
1645 : 0 : return
1646 : : !c
1647 : : !c last card of subroutine r1updt.
1648 : : !c
1649 : : end subroutine r1updt
1650 : :
1651 : : recursive &
1652 : 0 : &subroutine r1mpyq(m,n,a,lda,v,w)
1653 : : integer m,n,lda
1654 : : double precision a(lda,n),v(n),w(n)
1655 : : !c **********
1656 : : !c
1657 : : !c subroutine r1mpyq
1658 : : !c
1659 : : !c given an m by n matrix a, this subroutine computes a*q where
1660 : : !c q is the product of 2*(n - 1) transformations
1661 : : !c
1662 : : !c gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
1663 : : !c
1664 : : !c and gv(i), gw(i) are givens rotations in the (i,n) plane which
1665 : : !c eliminate elements in the i-th and n-th planes, respectively.
1666 : : !c q itself is not given, rather the information to recover the
1667 : : !c gv, gw rotations is supplied.
1668 : : !c
1669 : : !c the subroutine statement is
1670 : : !c
1671 : : !c subroutine r1mpyq(m,n,a,lda,v,w)
1672 : : !c
1673 : : !c where
1674 : : !c
1675 : : !c m is a positive integer input variable set to the number
1676 : : !c of rows of a.
1677 : : !c
1678 : : !c n is a positive integer input variable set to the number
1679 : : !c of columns of a.
1680 : : !c
1681 : : !c a is an m by n array. on input a must contain the matrix
1682 : : !c to be postmultiplied by the orthogonal matrix q
1683 : : !c described above. on output a*q has replaced a.
1684 : : !c
1685 : : !c lda is a positive integer input variable not less than m
1686 : : !c which specifies the leading dimension of the array a.
1687 : : !c
1688 : : !c v is an input array of length n. v(i) must contain the
1689 : : !c information necessary to recover the givens rotation gv(i)
1690 : : !c described above.
1691 : : !c
1692 : : !c w is an input array of length n. w(i) must contain the
1693 : : !c information necessary to recover the givens rotation gw(i)
1694 : : !c described above.
1695 : : !c
1696 : : !c subroutines called
1697 : : !c
1698 : : !c fortran-supplied ... dabs,dsqrt
1699 : : !c
1700 : : !c argonne national laboratory. minpack project. march 1980.
1701 : : !c burton s. garbow, kenneth e. hillstrom, jorge j. more
1702 : : !c
1703 : : !c **********
1704 : : integer i,j,nmj,nm1
1705 : : double precision cos,one,sin,temp
1706 : : data one /1.0d0/
1707 : : !c
1708 : : !c apply the first set of givens rotations to a.
1709 : : !c
1710 : 0 : nm1 = n - 1
1711 [ # # ]: 0 : if (nm1 .lt. 1) go to 50
1712 [ # # ]: 0 : do 20 nmj = 1, nm1
1713 : 0 : j = n - nmj
1714 [ # # ]: 0 : if (dabs(v(j)) .gt. one) cos = one/v(j)
1715 [ # # ]: 0 : if (dabs(v(j)) .gt. one) sin = dsqrt(one-cos**2)
1716 [ # # ]: 0 : if (dabs(v(j)) .le. one) sin = v(j)
1717 [ # # ]: 0 : if (dabs(v(j)) .le. one) cos = dsqrt(one-sin**2)
1718 [ # # ]: 0 : do 10 i = 1, m
1719 : 0 : temp = cos*a(i,j) - sin*a(i,n)
1720 : 0 : a(i,n) = sin*a(i,j) + cos*a(i,n)
1721 : 0 : a(i,j) = temp
1722 : 0 : 10 continue
1723 : 0 : 20 continue
1724 : : !c
1725 : : !c apply the second set of givens rotations to a.
1726 : : !c
1727 [ # # ]: 0 : do 40 j = 1, nm1
1728 [ # # ]: 0 : if (dabs(w(j)) .gt. one) cos = one/w(j)
1729 [ # # ]: 0 : if (dabs(w(j)) .gt. one) sin = dsqrt(one-cos**2)
1730 [ # # ]: 0 : if (dabs(w(j)) .le. one) sin = w(j)
1731 [ # # ]: 0 : if (dabs(w(j)) .le. one) cos = dsqrt(one-sin**2)
1732 [ # # ]: 0 : do 30 i = 1, m
1733 : 0 : temp = cos*a(i,j) + sin*a(i,n)
1734 : 0 : a(i,n) = -sin*a(i,j) + cos*a(i,n)
1735 : 0 : a(i,j) = temp
1736 : 0 : 30 continue
1737 : 0 : 40 continue
1738 : : 50 continue
1739 : 0 : return
1740 : : !c
1741 : : !c last card of subroutine r1mpyq.
1742 : : !c
1743 : : end subroutine r1mpyq
1744 : :
1745 : : !!***
1746 : :
1747 : : !----------------------------------------------------------------------
1748 : :
1749 : : END MODULE m_hybrd
1750 : : !!***
|
