← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/arpack/dgetv0.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/arpack/dgetv0.f
 
1
c-----------------------------------------------------------------------
 
2
c\BeginDoc
 
3
c
 
4
c\Name: dgetv0
 
5
c
 
6
c\Description: 
 
7
c  Generate a random initial residual vector for the Arnoldi process.
 
8
c  Force the residual vector to be in the range of the operator OP.  
 
9
c
 
10
c\Usage:
 
11
c  call dgetv0
 
12
c     ( IDO, BMAT, ITRY, INITV, N, J, V, LDV, RESID, RNORM, 
 
13
c       IPNTR, WORKD, IERR )
 
14
c
 
15
c\Arguments
 
16
c  IDO     Integer.  (INPUT/OUTPUT)
 
17
c          Reverse communication flag.  IDO must be zero on the first
 
18
c          call to dgetv0.
 
19
c          -------------------------------------------------------------
 
20
c          IDO =  0: first call to the reverse communication interface
 
21
c          IDO = -1: compute  Y = OP * X  where
 
22
c                    IPNTR(1) is the pointer into WORKD for X,
 
23
c                    IPNTR(2) is the pointer into WORKD for Y.
 
24
c                    This is for the initialization phase to force the
 
25
c                    starting vector into the range of OP.
 
26
c          IDO =  2: compute  Y = B * X  where
 
27
c                    IPNTR(1) is the pointer into WORKD for X,
 
28
c                    IPNTR(2) is the pointer into WORKD for Y.
 
29
c          IDO = 99: done
 
30
c          -------------------------------------------------------------
 
31
c
 
32
c  BMAT    Character*1.  (INPUT)
 
33
c          BMAT specifies the type of the matrix B in the (generalized)
 
34
c          eigenvalue problem A*x = lambda*B*x.
 
35
c          B = 'I' -> standard eigenvalue problem A*x = lambda*x
 
36
c          B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
 
37
c
 
38
c  ITRY    Integer.  (INPUT)
 
39
c          ITRY counts the number of times that dgetv0 is called.  
 
40
c          It should be set to 1 on the initial call to dgetv0.
 
41
c
 
42
c  INITV   Logical variable.  (INPUT)
 
43
c          .TRUE.  => the initial residual vector is given in RESID.
 
44
c          .FALSE. => generate a random initial residual vector.
 
45
c
 
46
c  N       Integer.  (INPUT)
 
47
c          Dimension of the problem.
 
48
c
 
49
c  J       Integer.  (INPUT)
 
50
c          Index of the residual vector to be generated, with respect to
 
51
c          the Arnoldi process.  J > 1 in case of a "restart".
 
52
c
 
53
c  V       Double precision N by J array.  (INPUT)
 
54
c          The first J-1 columns of V contain the current Arnoldi basis
 
55
c          if this is a "restart".
 
56
c
 
57
c  LDV     Integer.  (INPUT)
 
58
c          Leading dimension of V exactly as declared in the calling 
 
59
c          program.
 
60
c
 
61
c  RESID   Double precision array of length N.  (INPUT/OUTPUT)
 
62
c          Initial residual vector to be generated.  If RESID is 
 
63
c          provided, force RESID into the range of the operator OP.
 
64
c
 
65
c  RNORM   Double precision scalar.  (OUTPUT)
 
66
c          B-norm of the generated residual.
 
67
c
 
68
c  IPNTR   Integer array of length 3.  (OUTPUT)
 
69
c
 
70
c  WORKD   Double precision work array of length 2*N.  (REVERSE COMMUNICATION).
 
71
c          On exit, WORK(1:N) = B*RESID to be used in SSAITR.
 
72
c
 
73
c  IERR    Integer.  (OUTPUT)
 
74
c          =  0: Normal exit.
 
75
c          = -1: Cannot generate a nontrivial restarted residual vector
 
76
c                in the range of the operator OP.
 
77
c
 
78
c\EndDoc
 
79
c
 
80
c-----------------------------------------------------------------------
 
81
c
 
82
c\BeginLib
 
83
c
 
84
c\Local variables:
 
85
c     xxxxxx  real
 
86
c
 
87
c\References:
 
88
c  1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
 
89
c     a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
 
90
c     pp 357-385.
 
91
c  2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly 
 
92
c     Restarted Arnoldi Iteration", Rice University Technical Report
 
93
c     TR95-13, Department of Computational and Applied Mathematics.
 
94
c
 
95
c\Routines called:
 
96
c     second  ARPACK utility routine for timing.
 
97
c     dvout   ARPACK utility routine for vector output.
 
98
c     dlarnv  LAPACK routine for generating a random vector.
 
99
c     dgemv   Level 2 BLAS routine for matrix vector multiplication.
 
100
c     dcopy   Level 1 BLAS that copies one vector to another.
 
101
c     ddot    Level 1 BLAS that computes the scalar product of two vectors. 
 
102
c     dnrm2   Level 1 BLAS that computes the norm of a vector.
 
103
c
 
104
c\Author
 
105
c     Danny Sorensen               Phuong Vu
 
106
c     Richard Lehoucq              CRPC / Rice University
 
107
c     Dept. of Computational &     Houston, Texas
 
108
c     Applied Mathematics
 
109
c     Rice University           
 
110
c     Houston, Texas            
 
111
c
 
112
c\SCCS Information: @(#) 
 
113
c FILE: getv0.F   SID: 2.7   DATE OF SID: 04/07/99   RELEASE: 2
 
114
c
 
115
c\EndLib
 
116
c
 
117
c-----------------------------------------------------------------------
 
118
c
 
119
      subroutine dgetv0 
 
120
     &   ( ido, bmat, itry, initv, n, j, v, ldv, resid, rnorm, 
 
121
     &     ipntr, workd, ierr )
 
122
c 
 
123
c     %----------------------------------------------------%
 
124
c     | Include files for debugging and timing information |
 
125
c     %----------------------------------------------------%
 
126
c
 
127
      include   'debug.h'
 
128
      include   'stat.h'
 
129
c
 
130
c     %------------------%
 
131
c     | Scalar Arguments |
 
132
c     %------------------%
 
133
c
 
134
      character  bmat*1
 
135
      logical    initv
 
136
      integer    ido, ierr, itry, j, ldv, n
 
137
      Double precision
 
138
     &           rnorm
 
139
c
 
140
c     %-----------------%
 
141
c     | Array Arguments |
 
142
c     %-----------------%
 
143
c
 
144
      integer    ipntr(*)
 
145
      Double precision
 
146
     &           resid(n), v(ldv,j), workd(2*n)
 
147
c
 
148
c     %------------%
 
149
c     | Parameters |
 
150
c     %------------%
 
151
c
 
152
      Double precision
 
153
     &           one, zero
 
154
      parameter (one = 1.0D+0, zero = 0.0D+0)
 
155
c
 
156
c     %------------------------%
 
157
c     | Local Scalars & Arrays |
 
158
c     %------------------------%
 
159
c
 
160
      logical    first, inits, orth
 
161
      integer    idist, iseed(4), iter, msglvl, jj
 
162
      Double precision
 
163
     &           rnorm0
 
164
      save       first, iseed, inits, iter, msglvl, orth, rnorm0
 
165
c
 
166
c     %----------------------%
 
167
c     | External Subroutines |
 
168
c     %----------------------%
 
169
c
 
170
      external   dlarnv, dvout, dcopy, dgemv, second
 
171
c
 
172
c     %--------------------%
 
173
c     | External Functions |
 
174
c     %--------------------%
 
175
c
 
176
      Double precision
 
177
     &           ddot, dnrm2
 
178
      external   ddot, dnrm2
 
179
c
 
180
c     %---------------------%
 
181
c     | Intrinsic Functions |
 
182
c     %---------------------%
 
183
c
 
184
      intrinsic    abs, sqrt
 
185
c
 
186
c     %-----------------%
 
187
c     | Data Statements |
 
188
c     %-----------------%
 
189
c
 
190
      data       inits /.true./
 
191
c
 
192
c     %-----------------------%
 
193
c     | Executable Statements |
 
194
c     %-----------------------%
 
195
c
 
196
c
 
197
c     %-----------------------------------%
 
198
c     | Initialize the seed of the LAPACK |
 
199
c     | random number generator           |
 
200
c     %-----------------------------------%
 
201
c
 
202
      if (inits) then
 
203
          iseed(1) = 1
 
204
          iseed(2) = 3
 
205
          iseed(3) = 5
 
206
          iseed(4) = 7
 
207
          inits = .false.
 
208
      end if
 
209
c
 
210
      if (ido .eq.  0) then
 
211
c 
 
212
c        %-------------------------------%
 
213
c        | Initialize timing statistics  |
 
214
c        | & message level for debugging |
 
215
c        %-------------------------------%
 
216
c
 
217
         call second (t0)
 
218
         msglvl = mgetv0
 
219
c 
 
220
         ierr   = 0
 
221
         iter   = 0
 
222
         first  = .FALSE.
 
223
         orth   = .FALSE.
 
224
c
 
225
c        %-----------------------------------------------------%
 
226
c        | Possibly generate a random starting vector in RESID |
 
227
c        | Use a LAPACK random number generator used by the    |
 
228
c        | matrix generation routines.                         |
 
229
c        |    idist = 1: uniform (0,1)  distribution;          |
 
230
c        |    idist = 2: uniform (-1,1) distribution;          |
 
231
c        |    idist = 3: normal  (0,1)  distribution;          |
 
232
c        %-----------------------------------------------------%
 
233
c
 
234
         if (.not.initv) then
 
235
            idist = 2
 
236
            call dlarnv (idist, iseed, n, resid)
 
237
         end if
 
238
c 
 
239
c        %----------------------------------------------------------%
 
240
c        | Force the starting vector into the range of OP to handle |
 
241
c        | the generalized problem when B is possibly (singular).   |
 
242
c        %----------------------------------------------------------%
 
243
c
 
244
         call second (t2)
 
245
         if (bmat .eq. 'G') then
 
246
            nopx = nopx + 1
 
247
            ipntr(1) = 1
 
248
            ipntr(2) = n + 1
 
249
            call dcopy (n, resid, 1, workd, 1)
 
250
            ido = -1
 
251
            go to 9000
 
252
         end if
 
253
      end if
 
254
c 
 
255
c     %-----------------------------------------%
 
256
c     | Back from computing OP*(initial-vector) |
 
257
c     %-----------------------------------------%
 
258
c
 
259
      if (first) go to 20
 
260
c
 
261
c     %-----------------------------------------------%
 
262
c     | Back from computing B*(orthogonalized-vector) |
 
263
c     %-----------------------------------------------%
 
264
c
 
265
      if (orth)  go to 40
 
266
c 
 
267
      if (bmat .eq. 'G') then
 
268
         call second (t3)
 
269
         tmvopx = tmvopx + (t3 - t2)
 
270
      end if
 
271
c 
 
272
c     %------------------------------------------------------%
 
273
c     | Starting vector is now in the range of OP; r = OP*r; |
 
274
c     | Compute B-norm of starting vector.                   |
 
275
c     %------------------------------------------------------%
 
276
c
 
277
      call second (t2)
 
278
      first = .TRUE.
 
279
      if (bmat .eq. 'G') then
 
280
         nbx = nbx + 1
 
281
         call dcopy (n, workd(n+1), 1, resid, 1)
 
282
         ipntr(1) = n + 1
 
283
         ipntr(2) = 1
 
284
         ido = 2
 
285
         go to 9000
 
286
      else if (bmat .eq. 'I') then
 
287
         call dcopy (n, resid, 1, workd, 1)
 
288
      end if
 
289
c 
 
290
   20 continue
 
291
c
 
292
      if (bmat .eq. 'G') then
 
293
         call second (t3)
 
294
         tmvbx = tmvbx + (t3 - t2)
 
295
      end if
 
296
c 
 
297
      first = .FALSE.
 
298
      if (bmat .eq. 'G') then
 
299
          rnorm0 = ddot (n, resid, 1, workd, 1)
 
300
          rnorm0 = sqrt(abs(rnorm0))
 
301
      else if (bmat .eq. 'I') then
 
302
           rnorm0 = dnrm2(n, resid, 1)
 
303
      end if
 
304
      rnorm  = rnorm0
 
305
c
 
306
c     %---------------------------------------------%
 
307
c     | Exit if this is the very first Arnoldi step |
 
308
c     %---------------------------------------------%
 
309
c
 
310
      if (j .eq. 1) go to 50
 
311
c 
 
312
c     %----------------------------------------------------------------
 
313
c     | Otherwise need to B-orthogonalize the starting vector against |
 
314
c     | the current Arnoldi basis using Gram-Schmidt with iter. ref.  |
 
315
c     | This is the case where an invariant subspace is encountered   |
 
316
c     | in the middle of the Arnoldi factorization.                   |
 
317
c     |                                                               |
 
318
c     |       s = V^{T}*B*r;   r = r - V*s;                           |
 
319
c     |                                                               |
 
320
c     | Stopping criteria used for iter. ref. is discussed in         |
 
321
c     | Parlett's book, page 107 and in Gragg & Reichel TOMS paper.   |
 
322
c     %---------------------------------------------------------------%
 
323
c
 
324
      orth = .TRUE.
 
325
   30 continue
 
326
c
 
327
      call dgemv ('T', n, j-1, one, v, ldv, workd, 1, 
 
328
     &            zero, workd(n+1), 1)
 
329
      call dgemv ('N', n, j-1, -one, v, ldv, workd(n+1), 1, 
 
330
     &            one, resid, 1)
 
331
c 
 
332
c     %----------------------------------------------------------%
 
333
c     | Compute the B-norm of the orthogonalized starting vector |
 
334
c     %----------------------------------------------------------%
 
335
c
 
336
      call second (t2)
 
337
      if (bmat .eq. 'G') then
 
338
         nbx = nbx + 1
 
339
         call dcopy (n, resid, 1, workd(n+1), 1)
 
340
         ipntr(1) = n + 1
 
341
         ipntr(2) = 1
 
342
         ido = 2
 
343
         go to 9000
 
344
      else if (bmat .eq. 'I') then
 
345
         call dcopy (n, resid, 1, workd, 1)
 
346
      end if
 
347
c 
 
348
   40 continue
 
349
c
 
350
      if (bmat .eq. 'G') then
 
351
         call second (t3)
 
352
         tmvbx = tmvbx + (t3 - t2)
 
353
      end if
 
354
c 
 
355
      if (bmat .eq. 'G') then
 
356
         rnorm = ddot (n, resid, 1, workd, 1)
 
357
         rnorm = sqrt(abs(rnorm))
 
358
      else if (bmat .eq. 'I') then
 
359
         rnorm = dnrm2(n, resid, 1)
 
360
      end if
 
361
c
 
362
c     %--------------------------------------%
 
363
c     | Check for further orthogonalization. |
 
364
c     %--------------------------------------%
 
365
c
 
366
      if (msglvl .gt. 2) then
 
367
          call dvout (logfil, 1, rnorm0, ndigit, 
 
368
     &                '_getv0: re-orthonalization ; rnorm0 is')
 
369
          call dvout (logfil, 1, rnorm, ndigit, 
 
370
     &                '_getv0: re-orthonalization ; rnorm is')
 
371
      end if
 
372
c
 
373
      if (rnorm .gt. 0.717*rnorm0) go to 50
 
374
c 
 
375
      iter = iter + 1
 
376
      if (iter .le. 5) then
 
377
c
 
378
c        %-----------------------------------%
 
379
c        | Perform iterative refinement step |
 
380
c        %-----------------------------------%
 
381
c
 
382
         rnorm0 = rnorm
 
383
         go to 30
 
384
      else
 
385
c
 
386
c        %------------------------------------%
 
387
c        | Iterative refinement step "failed" |
 
388
c        %------------------------------------%
 
389
c
 
390
         do 45 jj = 1, n
 
391
            resid(jj) = zero
 
392
   45    continue
 
393
         rnorm = zero
 
394
         ierr = -1
 
395
      end if
 
396
c 
 
397
   50 continue
 
398
c
 
399
      if (msglvl .gt. 0) then
 
400
         call dvout (logfil, 1, rnorm, ndigit,
 
401
     &        '_getv0: B-norm of initial / restarted starting vector')
 
402
      end if
 
403
      if (msglvl .gt. 3) then
 
404
         call dvout (logfil, n, resid, ndigit,
 
405
     &        '_getv0: initial / restarted starting vector')
 
406
      end if
 
407
      ido = 99
 
408
c 
 
409
      call second (t1)
 
410
      tgetv0 = tgetv0 + (t1 - t0)
 
411
c 
 
412
 9000 continue
 
413
      return
 
414
c
 
415
c     %---------------%
 
416
c     | End of dgetv0 |
 
417
c     %---------------%
 
418
c
 
419
      end