SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pdhseqrdriver.f
Go to the documentation of this file.
1***********************************************************************
2* Test program for ScaLAPACK-style routine PDHSEQR *
3***********************************************************************
4*
5* Contributor: Robert Granat and Meiyue Shao
6* This version is of Feb 2011.
7*
9*
10* Declarations
11*
12 IMPLICIT NONE
13* ...Parameters...
14 LOGICAL balance, comphess, compresi,
15 $ comporth
16 LOGICAL debug, prn, timesteps, barr,
17 $ uni_lapack
18 INTEGER slv_min, slv_max
19 parameter( debug = .false.,
20 $ prn = .false.,
21 $ timesteps = .true.,
22 $ comphess = .true.,
23 $ compresi = .true.,
24 $ comporth = .true.,
25 $ balance = .true.,
26 $ barr = .false.,
27* Solver: 1-PDLAQR1, 2-PDHSEQR.
28 $ slv_min = 2, slv_max = 2,
29 $ uni_lapack = .true. )
30 INTEGER n, nb, arsrc, acsrc
31 parameter(
32* Problem size.
33 $ n = 500, nb = 50,
34* What processor should hold the first element in A?
35 $ arsrc = 0, acsrc = 0 )
36 INTEGER block_cyclic_2d, csrc_, ctxt_, dlen_, dt_,
37 $ lld_, mb_, m_, nb_, n_, rsrc_
38 parameter( block_cyclic_2d = 1, dlen_ = 9, dt_ = 1,
39 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
40 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
41 INTEGER dpalloc, intallc
42 INTEGER dpsiz, intsz, nout, izero
43 parameter( dpsiz = 8, dpalloc = 8 000 000,
44 $ intsz = 4, intallc = 8 000 000,
45 $ nout = 6, izero = 0)
46 DOUBLE PRECISION zero, one, two
47 parameter( zero = 0.0d+00, one = 1.0d+00, two = 2.0d+00 )
48*
49* ...Local Scalars...
50 INTEGER ictxt, iam, nprow, npcol, myrow, mycol,
51 $ sys_nprocs, nprocs, arows, acols, temp_ictxt
52 INTEGER threads
53 INTEGER info, ktop, kbot, ilo, ihi, i
54 INTEGER ipa, ipacpy, ipq, wr1, wi1, wr2, wi2, ipw1,
55 $ ipw2, ipiw
56 INTEGER totit, sweeps, totns, hess
57 DOUBLE PRECISION eps, thresh
58 DOUBLE PRECISION stamp, tottime, t_ba, t_gen, t_hs, t_sch, t_qr,
59 $ t_res, itpereig, swpspeig, nspeig, speedup,
60 $ efficiency
61 DOUBLE PRECISION rnorm, anorm, r1, orth, o1, o2, dpdum, elem1,
62 $ elem2, elem3, ediff
63 INTEGER solver
64 CHARACTER*6 passed
65*
66* ...Local Arrays...
67 INTEGER desca( dlen_ ), descq( dlen_ ), descvec( dlen_ )
68 DOUBLE PRECISION scale( n )
69 DOUBLE PRECISION, ALLOCATABLE :: mem(:)
70 INTEGER, ALLOCATABLE :: imem(:)
71*
72* ...Intrinsic Functions...
73 INTRINSIC int, dble, sqrt, max, min
74*
75* ...External Functions...
76 INTEGER numroc
77 DOUBLE PRECISION pdlamch, pdlange, mpi_wtime
78 EXTERNAL blacs_pinfo, blacs_get, blacs_gridinit,
79 $ blacs_gridinfo, blacs_gridexit, blacs_exit
81 EXTERNAL dgebal, dgehrd
82 EXTERNAL mpi_wtime
83 EXTERNAL pdgebal
84 EXTERNAL pdmatgen2
85*
86* ...Executable statements...
87*
88 CALL blacs_pinfo( iam, sys_nprocs )
89 nprow = int( sqrt( dble(sys_nprocs) ) )
90 npcol = sys_nprocs / nprow
91 CALL blacs_get( 0, 0, ictxt )
92 CALL blacs_gridinit( ictxt, '2D', nprow, npcol )
93 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
94c print*, iam, ictxt, myrow, mycol
95c IF ( ( MYROW.GE.NPROW ) .OR. ( MYCOL.GE.NPCOL ) ) GO TO 777
96 IF ( ictxt.LT.0 ) GO TO 777
97*
98* Read out the number of underlying threads and set stack size in
99* kilobytes.
100*
101 thresh = 30.0
102 tottime = mpi_wtime()
103 t_gen = 0.0d+00
104 t_res = 0.0d+00
105 t_sch = 0.0d+00
106*
107* Allocate and Init memory with zeros.
108*
109 info = 0
110 ALLOCATE ( mem( dpalloc ), stat = info )
111 IF( info.NE.0 ) THEN
112 WRITE(*,*) '% Could not allocate MEM. INFO = ', info
113 GO TO 777
114 END IF
115 ALLOCATE ( imem( intallc ), stat = info )
116 IF( info.NE.0 ) THEN
117 WRITE(*,*) '% Could not allocate IMEM. INFO = ', info
118 GO TO 777
119 END IF
120 mem( 1:dpalloc ) = zero
121 imem( 1:intallc ) = izero
122*
123* Get machine epsilon.
124*
125 eps = pdlamch( ictxt, 'Epsilon' )
126*
127* Print welcoming message.
128*
129 IF( iam.EQ.0 ) THEN
130 WRITE(*,*)
131 WRITE(*,*) 'ScaLAPACK Test for PDHSEQR'
132 WRITE(*,*)
133 WRITE(*,*) 'epsilon = ', eps
134 WRITE(*,*) 'threshold = ', thresh
135 WRITE(*,*)
136 WRITE(*,*) 'Residual and Orthogonality Residual computed by:'
137 WRITE(*,*)
138 WRITE(*,*) 'Residual = ',
139 $ ' || T - Q^T*A*Q ||_F / ( ||A||_F * eps * sqrt(N) )'
140 WRITE(*,*)
141 WRITE(*,*) 'Orthogonality = ',
142 $ ' MAX( || I - Q^T*Q ||_F, || I - Q*Q^T ||_F ) / ',
143 $ ' (eps * N)'
144 WRITE(*,*)
145 WRITE(*,*)
146 $ 'Test passes if both residuals are less then threshold'
147 WRITE( nout, * )
148 WRITE( nout, fmt = 9995 )
149 WRITE( nout, fmt = 9994 )
150 END IF
151*
152* Loop over problem parameters.
153*
154 DO ktop = 1, 1
155 DO kbot = n, n
156 DO solver = slv_max, slv_min, -1
157*
158* Set INFO to zero for this run.
159*
160 info = 0
161 nprocs = nprow*npcol
162 temp_ictxt = ictxt
163*
164* Count the number of rows and columns of current problem
165* for the current block sizes and grid properties.
166*
167 stamp = mpi_wtime()
168 arows = numroc( n, nb, myrow, 0, nprow )
169 acols = numroc( n, nb, mycol, 0, npcol )
170*
171* Set up matrix descriptors.
172*
173 IF( debug ) WRITE(*,*) '% #', iam, ': Set up descriptors...'
174 IF( barr ) CALL blacs_barrier(ictxt, 'A')
175 CALL descinit( desca, n, n, nb, nb, min(arsrc,nprow-1),
176 $ min(npcol-1,acsrc), temp_ictxt, max(1, arows), info )
177 IF ( info.NE.0 ) THEN
178 WRITE(*,*) "% DESCINIT DESCA failed, INFO =", info
179 GO TO 999
180 END IF
181 CALL descinit( descq, n, n, nb, nb, min(arsrc,nprow-1),
182 $ min(npcol-1,acsrc), temp_ictxt, max(1, arows), info )
183 IF ( info.NE.0 ) THEN
184 WRITE(*,*) "% DESCINIT DESCQ failed, INFO =", info
185 GO TO 999
186 END IF
187 CALL descinit( descvec, n, 1, n, 1, min(arsrc,nprow-1),
188 $ min(npcol-1,acsrc), temp_ictxt, n, info )
189 IF ( info.NE.0 ) THEN
190 WRITE(*,*) "% DESCINIT DESCVEC failed, INFO =", info
191 GO TO 999
192 END IF
193*
194* Assign pointer for ScaLAPACK arrays - first set DP memory.
195*
196 IF( debug ) WRITE(*,*) '% #', iam, ': Assign pointers...'
197 ipa = 1
198 ipacpy = ipa + desca( lld_ ) * acols
199 ipq = ipacpy + desca( lld_ ) * acols
200 wr1 = ipq + descq( lld_ ) * acols
201 wi1 = wr1 + n
202 wr2 = wi1 + n
203 wi2 = wr2 + n
204 ipw1 = wi2 + n
205 ipw2 = ipw1 + desca( lld_ ) * acols
206 IF( debug ) WRITE(*,*) '% (IPW2,DPALLOC):', ipw2, dpalloc
207* PRINT*, '%', IPA, IPACPY, IPQ, WR1, WI1, WR2, WI2, IPW1, IPW2
208 IF( ipw2+desca(lld_)*acols .GT. dpalloc+1 ) THEN
209 WRITE(*,*) '% Not enough DP memory!'
210 GO TO 999
211 END IF
212*
213* Then set integer memory pointers.
214*
215 ipiw = 1
216*
217* Generate testproblem.
218*
219 IF( barr ) CALL blacs_barrier(ictxt, 'A')
220 CALL pdlaset( 'All over', n, n, zero, one, mem(ipq), 1, 1,
221 $ descq )
222 CALL pdmatgen2( temp_ictxt, 'Random', 'NoDiagDominant',
223 $ n, n, nb, nb, mem(ipa), desca( lld_ ), 0, 0, 7, 0,
224 $ arows, 0, acols, myrow, mycol, nprow, npcol )
225 IF( .NOT. comphess ) THEN
226 CALL pdlaset( 'Lower triangular', n-2, n-2, zero, zero,
227 $ mem(ipa), 3, 1, desca )
228 CALL pdlaset( 'All over', n, n, zero, one, mem(ipq),
229 $ 1, 1, descq )
230 IF( ktop.GT.1 )
231 $ CALL pdlaset( 'Lower triangular', ktop-1, ktop-1,
232 $ zero, zero, mem(ipa), 2, 1, descq )
233 IF( kbot.LT.n )
234 $ CALL pdlaset( 'Lower triangular', n-kbot, n-kbot,
235 $ zero, zero, mem(ipa), kbot+1, kbot, descq )
236 END IF
237*
238* Do balancing if general matrix.
239*
240 IF( barr ) CALL blacs_barrier(ictxt, 'A')
241 t_ba = mpi_wtime()
242 IF( comphess .AND. balance ) THEN
243 IF( nprocs.EQ.1 .AND. solver.NE.2 .AND. uni_lapack ) THEN
244 IF( debug ) WRITE(*,*) '% #', iam, ': == dgebal =='
245 CALL dgebal( 'Both', n, mem(ipa), desca(lld_), ilo,
246 $ ihi, scale, info )
247 IF ( info.NE.0 ) THEN
248 WRITE(*,*) "% DGEBAL failed, INFO =", info
249 GO TO 999
250 END IF
251 ELSE
252 IF( debug ) WRITE(*,*) '% #', iam, ': == pdgebal =='
253 CALL pdgebal( 'Both', n, mem(ipa), desca, ilo, ihi,
254 $ scale, info )
255 IF ( info.NE.0 ) THEN
256 WRITE(*,*) "% PDGEBAL failed, INFO =", info
257 GO TO 999
258 END IF
259 END IF
260 ELSEIF( comphess ) THEN
261 ilo = 1
262 ihi = n
263 ELSE
264 ilo = ktop
265 ihi = kbot
266 END IF
267 t_ba = mpi_wtime() - t_ba
268c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
269c $ ' %%% Balancing took in seconds:',T_BA
270 IF( debug ) WRITE(*,*) '% #', iam, ': ILO,IHI=',ilo,ihi
271*
272* Make a copy of A.
273*
274 IF( barr ) CALL blacs_barrier(ictxt, 'A')
275 IF( debug ) WRITE(*,*) '% #', iam, ': Copy matrix A'
276 CALL pdlacpy( 'All', n, n, mem(ipa), 1, 1, desca, mem(ipacpy),
277 $ 1, 1, desca )
278*
279* Print matrices to screen in debugging mode.
280*
281 IF( prn )
282 $ CALL pdlaprnt( n, n, mem(ipacpy), 1, 1, desca, 0, 0,
283 $ 'A', nout, mem(ipw1) )
284 t_gen = t_gen + mpi_wtime() - stamp - t_ba
285c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
286c $ ' %%% Generation took in seconds:',T_GEN
287*
288* Only compute the Hessenberg form if necessary.
289*
290 t_hs = mpi_wtime()
291 IF( .NOT. comphess ) GO TO 30
292*
293* Reduce A to Hessenberg form.
294*
295 IF( barr ) CALL blacs_barrier(ictxt, 'A')
296 IF( debug ) WRITE(*,*) '% #', iam,
297 $ ': Reduce to Hessenberg form...N=',n, ilo,ihi
298* PRINT*, '% PDGEHRD: IPW2,MEM(IPW2)', IPW2, MEM(IPW2)
299 IF( nprocs.EQ.1 .AND. solver.NE.2 .AND. uni_lapack ) THEN
300 IF( debug ) WRITE(*,*) '% #', iam, ': == dgehrd =='
301 CALL dgehrd( n, ilo, ihi, mem(ipa), desca(lld_),
302 $ mem(ipw1), mem(ipw2), -1, info )
303 IF (dpalloc-ipw2.LT.mem(ipw2)) THEN
304 WRITE(*,*) "% Not enough memory for DGEHRD"
305 GO TO 999
306 END IF
307 CALL dgehrd( n, ilo, ihi, mem(ipa), desca(lld_),
308 $ mem(ipw1), mem(ipw2), dpalloc-ipw2, info )
309 IF ( info.NE.0 ) THEN
310 WRITE(*,*) "% DGEHRD failed, INFO =", info
311 GO TO 999
312 END IF
313 ELSE
314 IF( debug ) WRITE(*,*) '% #', iam, ': == pdgehrd =='
315 CALL pdgehrd( n, ilo, ihi, mem(ipa), 1, 1, desca, mem(ipw1),
316 $ mem(ipw2), -1, info )
317 IF (dpalloc-ipw2.LT.mem(ipw2)) THEN
318 WRITE(*,*) "% Not enough memory for PDGEHRD"
319 GO TO 999
320 END IF
321 CALL pdgehrd( n, ilo, ihi, mem(ipa), 1, 1, desca, mem(ipw1),
322 $ mem(ipw2), dpalloc-ipw2, info )
323 IF ( info.NE.0 ) THEN
324 WRITE(*,*) "% PDGEHRD failed, INFO =", info
325 GO TO 999
326 END IF
327 END IF
328*
329* Form Q explicitly.
330*
331 IF( barr ) CALL blacs_barrier(ictxt, 'A')
332 IF( debug ) WRITE(*,*) '% #', iam, ':Form Q explicitly'
333* PRINT*, '% PDORMHR: IPW2,MEM(IPW2)', IPW2, MEM(IPW2)
334 IF( debug ) WRITE(*,*) '% #', iam, ': == pdormhr =='
335 CALL pdormhr( 'L', 'N', n, n, ilo, ihi, mem(ipa), 1, 1,
336 $ desca, mem(ipw1), mem(ipq), 1, 1, descq, mem(ipw2),
337 $ -1, info )
338 IF (dpalloc-ipw2.LT.mem(ipw2)) THEN
339 WRITE(*,*) "% Not enough memory for PDORMHR"
340 GO TO 999
341 END IF
342 CALL pdormhr( 'L', 'N', n, n, ilo, ihi, mem(ipa), 1, 1,
343 $ desca, mem(ipw1), mem(ipq), 1, 1, descq, mem(ipw2),
344 $ dpalloc-ipw2, info )
345 IF ( info.NE.0 ) THEN
346 WRITE(*,*) "% PDORMHR failed, INFO =", info
347 GO TO 999
348 END IF
349*
350* Extract the upper Hessenberg part of A.
351*
352 CALL pdlaset( 'Lower triangular', n-2, n-2, zero, zero,
353 $ mem(ipa), 3, 1, desca )
354*
355* Print reduced matrix A in debugging mode.
356*
357 IF( prn ) THEN
358 CALL pdlaprnt( n, n, mem(ipa), 1, 1, desca, 0, 0, 'H', nout,
359 $ mem(ipw1) )
360 CALL pdlaprnt( n, n, mem(ipq), 1, 1, descq, 0, 0, 'Q', nout,
361 $ mem(ipw1) )
362 END IF
363*
364 30 CONTINUE
365 t_hs = mpi_wtime() - t_hs
366c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
367c $ ' %%% Hessenberg took in seconds:',T_HS
368*
369* Compute the real Schur form of the Hessenberg matrix A.
370*
371 IF( barr ) CALL blacs_barrier(ictxt, 'A')
372 t_qr = mpi_wtime()
373 IF( solver.EQ.1 ) THEN
374 IF( debug ) WRITE(*,*) '% #', iam, ': == pdlaqr1 =='
375* PRINT*, '% PDLAQR1: IPW1,MEM(IPW1)', IPW1, MEM(IPW1)
376 CALL pdlaqr1( .true., .true., n, ilo, ihi, mem(ipa), desca,
377 $ mem(wr1), mem(wi1), ilo, ihi, mem(ipq), descq,
378 $ mem(ipw1), -1, imem, -1, info )
379 IF (dpalloc-ipw1.LT.mem(ipw1)) THEN
380 WRITE(*,*) "% Not enough DP memory for PDLAQR1"
381 GO TO 999
382 END IF
383 IF (intallc.LT.imem(1)) THEN
384 WRITE(*,*) "% Not enough INT memory for PDLAQR1"
385 GO TO 999
386 END IF
387 CALL pdlaqr1( .true., .true., n, ilo, ihi, mem(ipa), desca,
388 $ mem(wr1), mem(wi1), ilo, ihi, mem(ipq), descq,
389 $ mem(ipw1), dpalloc-ipw1+1, imem, intallc, info )
390 IF (info.NE.0) THEN
391 WRITE(*,*) "% PDLAQR1: INFO =", info
392 END IF
393 ELSEIF( solver.EQ.2 ) THEN
394 IF( debug ) WRITE(*,*) '% #', iam, ': == pdhseqr =='
395* PRINT*, '% PDHSEQR: IPW1,MEM(IPW1)', IPW1, MEM(IPW1)
396 IF( barr ) CALL blacs_barrier(ictxt, 'A')
397 CALL pdhseqr( 'Schur', 'Vectors', n, ilo, ihi, mem(ipa),
398 $ desca, mem(wr2), mem(wi2), mem(ipq), descq, mem(ipw1),
399 $ -1, imem, -1, info )
400 IF( barr ) CALL blacs_barrier(ictxt, 'A')
401 IF (dpalloc-ipw1.LT.mem(ipw1)) THEN
402 WRITE(*,*) "% Not enough DP memory for PDHSEQR"
403 GO TO 999
404 END IF
405 IF (intallc.LT.imem(1)) THEN
406 WRITE(*,*) "% Not enough INT memory for PDHSEQR"
407 GO TO 999
408 END IF
409 IF( barr ) CALL blacs_barrier(ictxt, 'A')
410 CALL pdhseqr( 'Schur', 'Vectors', n, ilo, ihi, mem(ipa),
411 $ desca, mem(wr2), mem(wi2), mem(ipq), descq, mem(ipw1),
412 $ dpalloc-ipw1+1, imem, intallc, info )
413 IF (info.NE.0) THEN
414 WRITE(*,*) "% PDHSEQR: INFO =", info
415 END IF
416 ELSE
417 WRITE(*,*) '% ERROR: Illegal SOLVER number!'
418 GO TO 999
419 END IF
420 t_qr = mpi_wtime() - t_qr
421c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
422c $ ' %%% QR-algorithm took in seconds:',T_QR
423 t_sch = t_sch + t_qr + t_hs + t_ba
424* TOTIT = IMEM(1)
425* SWEEPS = IMEM(2)
426* TOTNS = IMEM(3)
427 itpereig = dble(totit) / dble(n)
428 swpspeig = dble(sweeps) / dble(n)
429 nspeig = dble(totns) / dble(n)
430*
431* Print reduced matrix A in debugging mode.
432*
433 IF( prn ) THEN
434 CALL pdlaprnt( n, n, mem(ipa), 1, 1, desca, 0, 0, 'T',
435 $ nout, mem(ipw1) )
436 CALL pdlaprnt( n, n, mem(ipq), 1, 1, descq, 0, 0, 'Z',
437 $ nout, mem(ipw1) )
438 END IF
439*
440* Check that returned Schur form is really a quasi-triangular
441* matrix.
442*
443 hess = 0
444 DO i = 1, n-1
445 IF( i.GT.1 ) THEN
446 CALL pdelget( 'All', '1-Tree', elem1, mem(ipa), i, i-1,
447 $ desca )
448 ELSE
449 elem1 = zero
450 END IF
451 CALL pdelget( 'All', '1-Tree', elem2, mem(ipa), i+1, i,
452 $ desca )
453 IF( i.LT.n-1 ) THEN
454 CALL pdelget( 'All', '1-Tree', elem3, mem(ipa), i+2, i+1,
455 $ desca )
456 ELSE
457 elem3 = zero
458 END IF
459 IF( elem2.NE.zero .AND. abs(elem1)+abs(elem2)+abs(elem3).GT.
460 $ abs(elem2) ) hess = hess + 1
461 END DO
462*
463* Compute residual norms and other results:
464*
465* 1) RNORM = || T - Q'*A*Q ||_F / ||A||_F
466* 2) ORTH = MAX( || I - Q'*Q ||_F, || I - Q*Q' ||_F ) /
467* (epsilon*N)
468*
469 stamp = mpi_wtime()
470 IF( compresi ) THEN
471 IF( debug ) WRITE(*,*) '% #', iam, ': Compute residuals 1'
472 IF( debug ) WRITE(*,*) '% #', iam, ': pdgemm 3'
473 CALL pdgemm( 'N', 'N', n, n, n, one, mem(ipacpy), 1, 1,
474 $ desca, mem(ipq), 1, 1, descq, zero, mem(ipw1), 1, 1,
475 $ desca )
476 IF( debug ) WRITE(*,*) '% #', iam, ': pdgemm 4'
477 IF( debug ) WRITE(*,*) '% #', iam, ': N=',n
478 IF( debug ) WRITE(*,*) '% #', iam, ': DESCA=',desca(1:dlen_)
479 IF( debug ) WRITE(*,*) '% #', iam, ': DESCQ=',descq(1:dlen_)
480 CALL pdgemm( 'T', 'N', n, n, n, -one, mem(ipq), 1, 1,
481 $ descq, mem(ipw1), 1, 1, desca, one, mem(ipa), 1, 1,
482 $ desca )
483 r1 = pdlange( 'Frobenius', n, n, mem(ipa), 1, 1, desca,
484 $ dpdum )
485 anorm = pdlange( 'Frobenius', n, n, mem(ipacpy), 1, 1,
486 $ desca, dpdum )
487 IF( anorm.GT.zero )THEN
488 rnorm = r1 / (anorm*eps*sqrt(dble(n)))
489 ELSE
490 rnorm = r1
491 END IF
492 ELSE
493 rnorm = 0.0d0
494 END IF
495*
496 IF( comporth ) THEN
497 IF( debug ) WRITE(*,*) '% #', iam, ': Compute residuals 2'
498 CALL pdlaset( 'All', n, n, zero, one, mem(ipw1), 1, 1,
499 $ descq )
500 CALL pdlacpy( 'All', n, n, mem(ipq), 1, 1, descq, mem(ipw2),
501 $ 1, 1, descq )
502 CALL pdgemm( 'T', 'N', n, n, n, -one, mem(ipq), 1, 1, descq,
503 $ mem(ipw2), 1, 1, descq, one, mem(ipw1), 1, 1, descq )
504 o1 = pdlange( 'Frobenius', n, n, mem(ipw1), 1, 1, descq,
505 $ dpdum )
506 CALL pdlaset( 'All', n, n, zero, one, mem(ipw1), 1, 1,
507 $ descq )
508 CALL pdgemm( 'N', 'T', n, n, n, -one, mem(ipq), 1, 1, descq,
509 $ mem(ipw2), 1, 1, descq, one, mem(ipw1), 1, 1, descq )
510 o2 = pdlange( 'Frobenius', n, n, mem(ipw1), 1, 1, descq,
511 $ dpdum )
512 orth = max(o1,o2) / (eps*dble(n))
513 ELSE
514 orth = 0.0d0
515 END IF
516*
517 t_res = t_res + mpi_wtime() - stamp
518c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
519c $ ' %%% Residuals took in seconds:',T_RES
520 tottime = mpi_wtime() - tottime
521c IF( IAM.EQ.0 ) WRITE(*,*)
522c $ ' %%% Total execution time in seconds:', TOTTIME
523*
524*
525* Print results to screen.
526*
527 IF( (orth.GT.thresh).OR.(rnorm.GT.thresh) ) THEN
528 passed = 'FAILED'
529 ELSE
530 passed = 'PASSED'
531 END IF
532 IF( debug ) WRITE(*,*) '% #', iam, ': Print results...'
533 IF( iam.EQ.0 ) THEN
534 WRITE( nout, fmt = 9993 ) n, nb, nprow, npcol, t_qr, passed
535 END IF
536 CALL blacs_barrier( ictxt, 'All' )
537 END DO
538 END DO
539 END DO
540 999 CONTINUE
541*
542* Deallocate MEM and IMEM.
543*
544 DEALLOCATE( mem, imem )
545*
546 CALL blacs_gridexit( ictxt )
547*
548 777 CONTINUE
549*
550 CALL blacs_exit( 0 )
551*
552* Format specifications.
553*
554 6666 FORMAT(a2,a3,a6,a4,a5,a6,a3,a3,a3,a9,a9,a9,a8,a8,a9,a9,a9,a9,a9,
555 $ a9,a9,a9,a9,a9,a9,a5,a5,a8,a5,a5)
556 7777 FORMAT(a2,i3,i6,i4,i5,i6,i3,i3,i3,f9.2,f9.2,f9.2,f8.2,f8.2,f9.2,
557 $ f9.2,f9.2,f9.2,f9.2,f9.2,f9.2,f9.2,e9.2,e9.2,e9.2,i5,i5,
558 $ f8.4,i5,i5,a2)
559 9995 FORMAT( ' N NB P Q QR Time CHECK' )
560 9994 FORMAT( '----- --- ---- ---- -------- ------' )
561 9993 FORMAT( i5, 1x, i3, 1x, i4, 1x, i4, 1x, f8.2, 1x, a6 )
562
563*
564 END
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
Definition descinit.f:3
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
#define max(A, B)
Definition pcgemr.c:180
#define min(A, B)
Definition pcgemr.c:181
subroutine pdlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pdblastst.f:6862
double precision function pdlamch(ictxt, cmach)
Definition pdblastst.f:6769
subroutine pdelget(scope, top, alpha, a, ia, ja, desca)
Definition pdelget.f:2
subroutine pdgebal(job, n, a, desca, ilo, ihi, scale, info)
Definition pdgebal.f:2
subroutine pdgehrd(n, ilo, ihi, a, ia, ja, desca, tau, work, lwork, info)
Definition pdgehrd.f:3
subroutine pdhseqr(job, compz, n, ilo, ihi, h, desch, wr, wi, z, descz, work, lwork, iwork, liwork, info)
Definition pdhseqr.f:3
program pdhseqrdriver
subroutine pdlacpy(uplo, m, n, a, ia, ja, desca, b, ib, jb, descb)
Definition pdlacpy.f:3
double precision function pdlange(norm, m, n, a, ia, ja, desca, work)
Definition pdlange.f:3
subroutine pdlaprnt(m, n, a, ia, ja, desca, irprnt, icprnt, cmatnm, nout, work)
Definition pdlaprnt.f:3
recursive subroutine pdlaqr1(wantt, wantz, n, ilo, ihi, a, desca, wr, wi, iloz, ihiz, z, descz, work, lwork, iwork, ilwork, info)
Definition pdlaqr1.f:5
subroutine pdmatgen2(ictxt, aform, diag, m, n, mb, nb, a, lda, iarow, iacol, iseed, iroff, irnum, icoff, icnum, myrow, mycol, nprow, npcol)
Definition pdmatgen2.f:4
subroutine pdormhr(side, trans, m, n, ilo, ihi, a, ia, ja, desca, tau, c, ic, jc, descc, work, lwork, info)
Definition pdormhr.f:3