1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
|
*********************************************************
*
* These routines are for testing LSQR.
*
*********************************************************
SUBROUTINE APROD ( MODE, M, N, X, Y, LENIW, LENRW, IW, RW )
INTEGER MODE, M, N, LENIW, LENRW
INTEGER IW(LENIW)
DOUBLE PRECISION X(N), Y(M), RW(LENRW)
* ------------------------------------------------------------------
* This is the matrix-vector product routine required by LSQR
* for a test matrix of the form A = HY*D*HZ. The quantities
* defining D, HY, HZ are in the work array RW, followed by a
* work array W. These are passed to APROD1 and APROD2 in order to
* make the code readable.
* ------------------------------------------------------------------
INTEGER LOCD, LOCHY, LOCHZ, LOCW
LOCD = 1
LOCHY = LOCD + N
LOCHZ = LOCHY + M
LOCW = LOCHZ + N
IF (MODE .EQ. 1) CALL APROD1( M, N, X, Y,
$ RW(LOCD), RW(LOCHY), RW(LOCHZ), RW(LOCW) )
IF (MODE .NE. 1) CALL APROD2( M, N, X, Y,
$ RW(LOCD), RW(LOCHY), RW(LOCHZ), RW(LOCW) )
* End of APROD
END
SUBROUTINE APROD1( M, N, X, Y, D, HY, HZ, W )
INTEGER M, N
DOUBLE PRECISION X(N), Y(M), D(N), HY(M), HZ(N), W(M)
* ------------------------------------------------------------------
* APROD1 computes Y = Y + A*X for subroutine APROD,
* where A is a test matrix of the form A = HY*D*HZ,
* and the latter matrices HY, D, HZ are represented by
* input vectors with the same name.
* ------------------------------------------------------------------
INTEGER I
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0 )
CALL HPROD ( N, HZ, X, W )
DO 100 I = 1, N
W(I) = D(I) * W(I)
100 CONTINUE
DO 200 I = N + 1, M
W(I) = ZERO
200 CONTINUE
CALL HPROD ( M, HY, W, W )
DO 600 I = 1, M
Y(I) = Y(I) + W(I)
600 CONTINUE
* End of APROD1
END
SUBROUTINE APROD2( M, N, X, Y, D, HY, HZ, W )
INTEGER M, N
DOUBLE PRECISION X(N), Y(M), D(N), HY(M), HZ(N), W(M)
* ------------------------------------------------------------------
* APROD2 computes X = X + A(T)*Y for subroutine APROD,
* where A is a test matrix of the form A = HY*D*HZ,
* and the latter matrices HY, D, HZ are represented by
* input vectors with the same name.
* ------------------------------------------------------------------
INTEGER I
CALL HPROD ( M, HY, Y, W )
DO 100 I = 1, N
W(I) = D(I)*W(I)
100 CONTINUE
CALL HPROD ( N, HZ, W, W )
DO 600 I = 1, N
X(I) = X(I) + W(I)
600 CONTINUE
* End of APROD2
END
SUBROUTINE HPROD ( N, HZ, X, Y )
INTEGER N
DOUBLE PRECISION HZ(N), X(N), Y(N)
* ------------------------------------------------------------------
* HPROD applies a Householder transformation stored in HZ
* to get Y = ( I - 2*HZ*HZ(transpose) ) * X.
* ------------------------------------------------------------------
INTEGER I
DOUBLE PRECISION S
S = 0.0
DO 100 I = 1, N
S = HZ(I) * X(I) + S
100 CONTINUE
S = S + S
DO 200 I = 1, N
Y(I) = X(I) - S * HZ(I)
200 CONTINUE
* End of HPROD
END
SUBROUTINE LSTP ( M, N, NDUPLC, NPOWER, DAMP, X,
$ B, D, HY, HZ, W, ACOND, RNORM )
INTEGER M, N, MAXMN, NDUPLC, NPOWER
DOUBLE PRECISION DAMP, ACOND, RNORM
DOUBLE PRECISION B(M), X(N), D(N), HY(M), HZ(N), W(M)
* ------------------------------------------------------------------
* LSTP generates a sparse least-squares test problem of the form
*
* ( A )*X = ( B )
* ( DAMP*I ) ( 0 )
*
* having a specified solution X. The matrix A is constructed
* in the form A = HY*D*HZ, where D is an M by N diagonal matrix,
* and HY and HZ are Householder transformations.
*
* The first 6 parameters are input to LSTP. The remainder are
* output. LSTP is intended for use with M .GE. N.
*
*
* Functions and subroutines
*
* TESTPROB APROD1, HPROD
* BLAS DNRM2
* ------------------------------------------------------------------
* Intrinsics and local variables
INTRINSIC COS, SIN, SQRT
INTEGER I, J
DOUBLE PRECISION DNRM2
DOUBLE PRECISION ALFA, BETA, DAMPSQ, FOURPI, T
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0, ONE = 1.0 )
* ------------------------------------------------------------------
* Make two vectors of norm 1.0 for the Householder transformations.
* FOURPI need not be exact.
* ------------------------------------------------------------------
DAMPSQ = DAMP**2
FOURPI = 4.0 * 3.141592
ALFA = FOURPI / M
BETA = FOURPI / N
DO 100 I = 1, M
HY(I) = SIN( I * ALFA )
100 CONTINUE
DO 200 I = 1, N
HZ(I) = COS( I * BETA )
200 CONTINUE
ALFA = DNRM2 ( M, HY, 1 )
BETA = DNRM2 ( N, HZ, 1 )
CALL DSCAL ( M, (- ONE / ALFA), HY, 1 )
CALL DSCAL ( N, (- ONE / BETA), HZ, 1 )
*
* ------------------------------------------------------------------
* Set the diagonal matrix D. These are the singular values of A.
* ------------------------------------------------------------------
DO 300 I = 1, N
J = (I - 1 + NDUPLC) / NDUPLC
T = J * NDUPLC
T = T / N
D(I) = T**NPOWER
300 CONTINUE
ACOND = SQRT( (D(N)**2 + DAMPSQ) / (D(1)**2 + DAMPSQ) )
* ------------------------------------------------------------------
* Compute the residual vector, storing it in B.
* It takes the form HY*( s )
* ( t )
* where s is obtained from D*s = DAMP**2 * HZ * X
* and t can be anything.
* ------------------------------------------------------------------
CALL HPROD ( N, HZ, X, B )
DO 500 I = 1, N
B(I) = DAMPSQ * B(I) / D(I)
500 CONTINUE
T = ONE
DO 600 I = N + 1, M
J = I - N
B(I) = (T * J) / M
T = - T
600 CONTINUE
CALL HPROD ( M, HY, B, B )
* ------------------------------------------------------------------
* Now compute the true B = RESIDUAL + A*X.
* ------------------------------------------------------------------
RNORM = SQRT( DNRM2 ( M, B, 1 )**2
$ + DAMPSQ * DNRM2 ( N, X, 1 )**2 )
CALL APROD1( M, N, X, B, D, HY, HZ, W )
* End of LSTP
END
SUBROUTINE TEST ( M, N, NDUPLC, NPOWER, DAMP )
INTEGER M, N, NDUPLC, NPOWER
DOUBLE PRECISION DAMP
* ------------------------------------------------------------------
* This is an example driver routine for running LSQR.
* It generates a test problem, solves it, and examines the results.
* Note that subroutine APROD must be declared EXTERNAL
* if it is used only in the call to LSQR.
*
*
* Functions and subroutines
*
* TESTPROB APROD
* BLAS DCOPY, DNRM2, DSCAL
* ------------------------------------------------------------------
* Intrinsics and local variables
INTRINSIC MAX, SQRT
EXTERNAL APROD
INTEGER ISTOP, ITNLIM, J, NOUT
DOUBLE PRECISION DNRM2
PARAMETER ( MAXM = 200, MAXN = 100 )
DOUBLE PRECISION B(MAXM), U(MAXM),
$ V(MAXN), W(MAXN), X(MAXN),
$ SE(MAXN), XTRUE(MAXN)
DOUBLE PRECISION ATOL, BTOL, CONLIM,
$ ANORM, ACOND, RNORM, ARNORM,
$ DAMPSQ, ENORM, ETOL, XNORM
PARAMETER ( LENIW = 1, LENRW = 600 )
INTEGER IW(LENIW)
DOUBLE PRECISION RW(LENRW)
INTEGER LOCD, LOCHY, LOCHZ, LOCW, LTOTAL
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0 )
CHARACTER*34 LINE
DATA LINE
$ /'----------------------------------'/
* Set the desired solution XTRUE.
DO 100 J = 1, N
XTRUE(J) = N - J
100 CONTINUE
* Generate the specified test problem.
* The workspace array IW is not needed in this application.
* The workspace array RW is used for the following vectors:
* D(N), HY(M), HZ(N), W(MAX(M,N)).
* The vectors D, HY, HZ will define the test matrix A.
* W is needed for workspace in APROD1 and APROD2.
LOCD = 1
LOCHY = LOCD + N
LOCHZ = LOCHY + M
LOCW = LOCHZ + N
LTOTAL = LOCW + MAX(M,N) - 1
IF (LTOTAL .GT. LENRW) GO TO 900
CALL LSTP ( M, N, NDUPLC, NPOWER, DAMP, XTRUE,
$ B, RW(LOCD), RW(LOCHY), RW(LOCHZ), RW(LOCW),
$ ACOND, RNORM )
* Solve the problem defined by APROD, DAMP and B.
* Copy the rhs vector B into U (LSQR will overwrite U)
* and set the other input parameters for LSQR.
CALL DCOPY ( M, B, 1, U, 1 )
ATOL = 1.0E-10
BTOL = ATOL
CONLIM = 10.0 * ACOND
ITNLIM = M + N + 50
NOUT = 6
WRITE(NOUT, 1000) LINE, LINE,
$ M, N, NDUPLC, NPOWER, DAMP, ACOND, RNORM,
$ LINE, LINE
CALL LSQR ( M, N, APROD, DAMP,
$ LENIW, LENRW, IW, RW,
$ U, V, W, X, SE,
$ ATOL, BTOL, CONLIM, ITNLIM, NOUT,
$ ISTOP, ITN, ANORM, ACOND, RNORM, ARNORM, XNORM )
* Examine the results.
* We print the residual norms RNORM and ARNORM given by LSQR,
* and then compute their true values in terms of the solution X
* obtained by LSQR. At least one of them should be small.
DAMPSQ = DAMP**2
WRITE(NOUT, 2000)
WRITE(NOUT, 2100) RNORM, ARNORM, XNORM
* Compute U = A*X - B.
* This is the negative of the usual residual vector.
* It will be close to zero only if B is a compatible rhs
* and X is close to a solution.
CALL DCOPY ( M, B, 1, U, 1 )
CALL DSCAL ( M, (-ONE), U, 1 )
CALL APROD ( 1, M, N, X, U, LENIW, LENRW, IW, RW )
* Compute V = A(transpose)*U + DAMP**2 * X.
* This will be close to zero in all cases
* if X is close to a solution.
CALL DCOPY ( N, X, 1, V, 1 )
CALL DSCAL ( N, DAMPSQ, V, 1 )
CALL APROD ( 2, M, N, V, U, LENIW, LENRW, IW, RW )
* Compute the norms associated with X, U, V.
XNORM = DNRM2 ( N, X, 1 )
RNORM = SQRT( DNRM2 ( M, U, 1 )**2 + DAMPSQ * XNORM**2 )
ARNORM = DNRM2 ( N, V, 1 )
WRITE(NOUT, 2200) RNORM, ARNORM, XNORM
* Print the solution and standard error estimates from LSQR.
WRITE(NOUT, 2500) (J, X(J), J = 1, N)
WRITE(NOUT, 2600) (J, SE(J), J = 1, N)
* Print a clue about whether the solution looks OK.
DO 500 J = 1, N
W(J) = X(J) - XTRUE(J)
500 CONTINUE
ENORM = DNRM2 ( N, W, 1 ) / (ONE + DNRM2 ( N, XTRUE, 1 ))
ETOL = 1.0D-5
IF (ENORM .LE. ETOL) WRITE(NOUT, 3000) ENORM
IF (ENORM .GT. ETOL) WRITE(NOUT, 3100) ENORM
RETURN
* Not enough workspace.
900 WRITE(NOUT, 9000) LTOTAL
RETURN
1000 FORMAT(1P
$ // 1X, 2A
$ / ' Least-Squares Test Problem P(', 4I5, E12.2, ' )'
$ // ' Condition no. =', E12.4, ' Residual function =', E17.9
$ / 1X, 2A)
2000 FORMAT(
$ // 22X, ' Residual norm Residual norm Solution norm'
$ / 22X, '(Abar X - bbar) (Normal eqns) (X)' /)
2100 FORMAT(1P, ' Estimated by LSQR', 3E17.5)
2200 FORMAT(1P, ' Computed from X ', 3E17.5)
2500 FORMAT(//' Solution X' / 4(I6, G14.6))
2600 FORMAT(/ ' Standard errors SE' / 4(I6, G14.6))
3000 FORMAT(1P / ' LSQR appears to be successful.',
$ ' Relative error in X =', E10.2)
3100 FORMAT(1P / ' LSQR appears to have failed. ',
$ ' Relative error in X =', E10.2)
9000 FORMAT(/ ' XXX Insufficient workspace.',
$ ' The length of RW should be at least', I6)
* End of TEST
END
* -------------
* Main program.
* -------------
DOUBLE PRECISION DAMP1, DAMP2, DAMP3, DAMP4, ZERO
*
ZERO = 0.0
DAMP1 = 0.1
DAMP2 = 0.01
DAMP3 = 0.001
DAMP4 = 0.0001
CALL TEST ( 1, 1, 1, 1, ZERO )
CALL TEST ( 2, 1, 1, 1, ZERO )
CALL TEST ( 40, 40, 4, 4, ZERO )
CALL TEST ( 40, 40, 4, 4, DAMP2 )
CALL TEST ( 80, 40, 4, 4, DAMP2 )
STOP
* End of main program for testing LSQR
END
|
