ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
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 *
8  PROGRAM pdhseqrdriver
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
80  EXTERNAL numroc, pdlamch, pdlaset, pdgehrd, pdlange
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 )
94 c print*, iam, ictxt, myrow, mycol
95 c 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
268 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
269 c $ ' %%% 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
285 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
286 c $ ' %%% 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
366 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
367 c $ ' %%% 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
421 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
422 c $ ' %%% 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
518 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
519 c $ ' %%% Residuals took in seconds:',T_RES
520  tottime = mpi_wtime() - tottime
521 c IF( IAM.EQ.0 ) WRITE(*,*)
522 c $ ' %%% 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
max
#define max(A, B)
Definition: pcgemr.c:180
pdormhr
subroutine pdormhr(SIDE, TRANS, M, N, ILO, IHI, A, IA, JA, DESCA, TAU, C, IC, JC, DESCC, WORK, LWORK, INFO)
Definition: pdormhr.f:3
pdelget
subroutine pdelget(SCOPE, TOP, ALPHA, A, IA, JA, DESCA)
Definition: pdelget.f:2
pdlaqr1
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
pdgebal
subroutine pdgebal(JOB, N, A, DESCA, ILO, IHI, SCALE, INFO)
Definition: pdgebal.f:2
pdgehrd
subroutine pdgehrd(N, ILO, IHI, A, IA, JA, DESCA, TAU, WORK, LWORK, INFO)
Definition: pdgehrd.f:3
pdhseqr
subroutine pdhseqr(JOB, COMPZ, N, ILO, IHI, H, DESCH, WR, WI, Z, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdhseqr.f:3
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
pdlange
double precision function pdlange(NORM, M, N, A, IA, JA, DESCA, WORK)
Definition: pdlange.f:3
pdmatgen2
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
pdlaset
subroutine pdlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pdblastst.f:6862
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
pdlaprnt
subroutine pdlaprnt(M, N, A, IA, JA, DESCA, IRPRNT, ICPRNT, CMATNM, NOUT, WORK)
Definition: pdlaprnt.f:3
pdhseqrdriver
program pdhseqrdriver
Definition: pdhseqrdriver.f:8
pdlacpy
subroutine pdlacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pdlacpy.f:3
pdlamch
double precision function pdlamch(ICTXT, CMACH)
Definition: pdblastst.f:6769
min
#define min(A, B)
Definition: pcgemr.c:181