ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pshseqrdriver.f
Go to the documentation of this file.
1 ***********************************************************************
2 * Test program for ScaLAPACK-style routine PSHSEQR *
3 ***********************************************************************
4 *
5 * Contributor: Robert Granat and Meiyue Shao
6 * This version is of Feb 2011.
7 *
8  PROGRAM pshseqrdriver
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-PSLAQR1, 2-PSHSEQR.
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  REAL zero, one, two
47  parameter( zero = 0.0, one = 1.0, two = 2.0 )
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  REAL 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  REAL 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  REAL scale( n )
69  REAL, ALLOCATABLE :: mem(:)
70  INTEGER, ALLOCATABLE :: imem(:)
71 *
72 * ...Intrinsic Functions...
73  INTRINSIC int, float, sqrt, max, min
74 *
75 * ...External Functions...
76  INTEGER numroc
77  REAL pslamch, pslange
78  DOUBLE PRECISION mpi_wtime
79  EXTERNAL blacs_pinfo, blacs_get, blacs_gridinit,
80  $ blacs_gridinfo, blacs_gridexit, blacs_exit
81  EXTERNAL numroc, pslamch, pslaset, psgehrd, pslange
82  EXTERNAL sgebal, sgehrd
83  EXTERNAL mpi_wtime
84  EXTERNAL psgebal
85  EXTERNAL psmatgen2
86 *
87 * ...Executable statements...
88 *
89  CALL blacs_pinfo( iam, sys_nprocs )
90  nprow = int( sqrt( float(sys_nprocs) ) )
91  npcol = sys_nprocs / nprow
92  CALL blacs_get( 0, 0, ictxt )
93  CALL blacs_gridinit( ictxt, '2D', nprow, npcol )
94  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
95 c print*, iam, ictxt, myrow, mycol
96 c IF ( ( MYROW.GE.NPROW ) .OR. ( MYCOL.GE.NPCOL ) ) GO TO 777
97  IF ( ictxt.LT.0 ) GO TO 777
98 *
99 * Read out the number of underlying threads and set stack size in
100 * kilobytes.
101 *
102  thresh = 30.0
103  tottime = mpi_wtime()
104  t_gen = 0.0d+00
105  t_res = 0.0d+00
106  t_sch = 0.0d+00
107 *
108 * Allocate and Init memory with zeros.
109 *
110  info = 0
111  ALLOCATE ( mem( dpalloc ), stat = info )
112  IF( info.NE.0 ) THEN
113  WRITE(*,*) '% Could not allocate MEM. INFO = ', info
114  GO TO 777
115  END IF
116  ALLOCATE ( imem( intallc ), stat = info )
117  IF( info.NE.0 ) THEN
118  WRITE(*,*) '% Could not allocate IMEM. INFO = ', info
119  GO TO 777
120  END IF
121  mem( 1:dpalloc ) = zero
122  imem( 1:intallc ) = izero
123 *
124 * Get machine epsilon.
125 *
126  eps = pslamch( ictxt, 'Epsilon' )
127 *
128 * Print welcoming message.
129 *
130  IF( iam.EQ.0 ) THEN
131  WRITE(*,*)
132  WRITE(*,*) 'ScaLAPACK Test for PSHSEQR'
133  WRITE(*,*)
134  WRITE(*,*) 'epsilon = ', eps
135  WRITE(*,*) 'threshold = ', thresh
136  WRITE(*,*)
137  WRITE(*,*) 'Residual and Orthogonality Residual computed by:'
138  WRITE(*,*)
139  WRITE(*,*) 'Residual = ',
140  $ ' || T - Q^T*A*Q ||_F / ( ||A||_F * eps * sqrt(N) )'
141  WRITE(*,*)
142  WRITE(*,*) 'Orthogonality = ',
143  $ ' MAX( || I - Q^T*Q ||_F, || I - Q*Q^T ||_F ) / ',
144  $ ' (eps * N)'
145  WRITE(*,*)
146  WRITE(*,*)
147  $ 'Test passes if both residuals are less then threshold'
148  WRITE( nout, * )
149  WRITE( nout, fmt = 9995 )
150  WRITE( nout, fmt = 9994 )
151  END IF
152 *
153 * Loop over problem parameters.
154 *
155  DO ktop = 1, 1
156  DO kbot = n, n
157  DO solver = slv_max, slv_min, -1
158 *
159 * Set INFO to zero for this run.
160 *
161  info = 0
162  nprocs = nprow*npcol
163  temp_ictxt = ictxt
164 *
165 * Count the number of rows and columns of current problem
166 * for the current block sizes and grid properties.
167 *
168  stamp = mpi_wtime()
169  arows = numroc( n, nb, myrow, 0, nprow )
170  acols = numroc( n, nb, mycol, 0, npcol )
171 *
172 * Set up matrix descriptors.
173 *
174  IF( debug ) WRITE(*,*) '% #', iam, ': Set up descriptors...'
175  IF( barr ) CALL blacs_barrier(ictxt, 'A')
176  CALL descinit( desca, n, n, nb, nb, min(arsrc,nprow-1),
177  $ min(npcol-1,acsrc), temp_ictxt, max(1, arows), info )
178  IF ( info.NE.0 ) THEN
179  WRITE(*,*) "% DESCINIT DESCA failed, INFO =", info
180  GO TO 999
181  END IF
182  CALL descinit( descq, n, n, nb, nb, min(arsrc,nprow-1),
183  $ min(npcol-1,acsrc), temp_ictxt, max(1, arows), info )
184  IF ( info.NE.0 ) THEN
185  WRITE(*,*) "% DESCINIT DESCQ failed, INFO =", info
186  GO TO 999
187  END IF
188  CALL descinit( descvec, n, 1, n, 1, min(arsrc,nprow-1),
189  $ min(npcol-1,acsrc), temp_ictxt, n, info )
190  IF ( info.NE.0 ) THEN
191  WRITE(*,*) "% DESCINIT DESCVEC failed, INFO =", info
192  GO TO 999
193  END IF
194 *
195 * Assign pointer for ScaLAPACK arrays - first set DP memory.
196 *
197  IF( debug ) WRITE(*,*) '% #', iam, ': Assign pointers...'
198  ipa = 1
199  ipacpy = ipa + desca( lld_ ) * acols
200  ipq = ipacpy + desca( lld_ ) * acols
201  wr1 = ipq + descq( lld_ ) * acols
202  wi1 = wr1 + n
203  wr2 = wi1 + n
204  wi2 = wr2 + n
205  ipw1 = wi2 + n
206  ipw2 = ipw1 + desca( lld_ ) * acols
207  IF( debug ) WRITE(*,*) '% (IPW2,DPALLOC):', ipw2, dpalloc
208 * PRINT*, '%', IPA, IPACPY, IPQ, WR1, WI1, WR2, WI2, IPW1, IPW2
209  IF( ipw2+desca(lld_)*acols .GT. dpalloc+1 ) THEN
210  WRITE(*,*) '% Not enough DP memory!'
211  GO TO 999
212  END IF
213 *
214 * Then set integer memory pointers.
215 *
216  ipiw = 1
217 *
218 * Generate testproblem.
219 *
220  IF( barr ) CALL blacs_barrier(ictxt, 'A')
221  CALL pslaset( 'All over', n, n, zero, one, mem(ipq), 1, 1,
222  $ descq )
223  CALL psmatgen2( temp_ictxt, 'Random', 'NoDiagDominant',
224  $ n, n, nb, nb, mem(ipa), desca( lld_ ), 0, 0, 7, 0,
225  $ arows, 0, acols, myrow, mycol, nprow, npcol )
226  IF( .NOT. comphess ) THEN
227  CALL pslaset( 'Lower triangular', n-2, n-2, zero, zero,
228  $ mem(ipa), 3, 1, desca )
229  CALL pslaset( 'All over', n, n, zero, one, mem(ipq),
230  $ 1, 1, descq )
231  IF( ktop.GT.1 )
232  $ CALL pslaset( 'Lower triangular', ktop-1, ktop-1,
233  $ zero, zero, mem(ipa), 2, 1, descq )
234  IF( kbot.LT.n )
235  $ CALL pslaset( 'Lower triangular', n-kbot, n-kbot,
236  $ zero, zero, mem(ipa), kbot+1, kbot, descq )
237  END IF
238 *
239 * Do balancing if general matrix.
240 *
241  IF( barr ) CALL blacs_barrier(ictxt, 'A')
242  t_ba = mpi_wtime()
243  IF( comphess .AND. balance ) THEN
244  IF( nprocs.EQ.1 .AND. solver.NE.2 .AND. uni_lapack ) THEN
245  IF( debug ) WRITE(*,*) '% #', iam, ': == dgebal =='
246  CALL sgebal( 'Both', n, mem(ipa), desca(lld_), ilo,
247  $ ihi, scale, info )
248  IF ( info.NE.0 ) THEN
249  WRITE(*,*) "% SGEBAL failed, INFO =", info
250  GO TO 999
251  END IF
252  ELSE
253  IF( debug ) WRITE(*,*) '% #', iam, ': == pdgebal =='
254  CALL psgebal( 'Both', n, mem(ipa), desca, ilo, ihi,
255  $ scale, info )
256  IF ( info.NE.0 ) THEN
257  WRITE(*,*) "% PSGEBAL failed, INFO =", info
258  GO TO 999
259  END IF
260  END IF
261  ELSEIF( comphess ) THEN
262  ilo = 1
263  ihi = n
264  ELSE
265  ilo = ktop
266  ihi = kbot
267  END IF
268  t_ba = mpi_wtime() - t_ba
269 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
270 c $ ' %%% Balancing took in seconds:',T_BA
271  IF( debug ) WRITE(*,*) '% #', iam, ': ILO,IHI=',ilo,ihi
272 *
273 * Make a copy of A.
274 *
275  IF( barr ) CALL blacs_barrier(ictxt, 'A')
276  IF( debug ) WRITE(*,*) '% #', iam, ': Copy matrix A'
277  CALL pslacpy( 'All', n, n, mem(ipa), 1, 1, desca, mem(ipacpy),
278  $ 1, 1, desca )
279 *
280 * Print matrices to screen in debugging mode.
281 *
282  IF( prn )
283  $ CALL pslaprnt( n, n, mem(ipacpy), 1, 1, desca, 0, 0,
284  $ 'A', nout, mem(ipw1) )
285  t_gen = t_gen + mpi_wtime() - stamp - t_ba
286 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
287 c $ ' %%% Generation took in seconds:',T_GEN
288 *
289 * Only compute the Hessenberg form if necessary.
290 *
291  t_hs = mpi_wtime()
292  IF( .NOT. comphess ) GO TO 30
293 *
294 * Reduce A to Hessenberg form.
295 *
296  IF( barr ) CALL blacs_barrier(ictxt, 'A')
297  IF( debug ) WRITE(*,*) '% #', iam,
298  $ ': Reduce to Hessenberg form...N=',n, ilo,ihi
299 * PRINT*, '% PSGEHRD: IPW2,MEM(IPW2)', IPW2, MEM(IPW2)
300  IF( nprocs.EQ.1 .AND. solver.NE.2 .AND. uni_lapack ) THEN
301  IF( debug ) WRITE(*,*) '% #', iam, ': == dgehrd =='
302  CALL sgehrd( n, ilo, ihi, mem(ipa), desca(lld_),
303  $ mem(ipw1), mem(ipw2), -1, info )
304  IF (dpalloc-ipw2.LT.mem(ipw2)) THEN
305  WRITE(*,*) "% Not enough memory for SGEHRD"
306  GO TO 999
307  END IF
308  CALL sgehrd( n, ilo, ihi, mem(ipa), desca(lld_),
309  $ mem(ipw1), mem(ipw2), dpalloc-ipw2, info )
310  IF ( info.NE.0 ) THEN
311  WRITE(*,*) "% SGEHRD failed, INFO =", info
312  GO TO 999
313  END IF
314  ELSE
315  IF( debug ) WRITE(*,*) '% #', iam, ': == pdgehrd =='
316  CALL psgehrd( n, ilo, ihi, mem(ipa), 1, 1, desca, mem(ipw1),
317  $ mem(ipw2), -1, info )
318  IF (dpalloc-ipw2.LT.mem(ipw2)) THEN
319  WRITE(*,*) "% Not enough memory for PSGEHRD"
320  GO TO 999
321  END IF
322  CALL psgehrd( n, ilo, ihi, mem(ipa), 1, 1, desca, mem(ipw1),
323  $ mem(ipw2), dpalloc-ipw2, info )
324  IF ( info.NE.0 ) THEN
325  WRITE(*,*) "% PSGEHRD failed, INFO =", info
326  GO TO 999
327  END IF
328  END IF
329 *
330 * Form Q explicitly.
331 *
332  IF( barr ) CALL blacs_barrier(ictxt, 'A')
333  IF( debug ) WRITE(*,*) '% #', iam, ':Form Q explicitly'
334 * PRINT*, '% PSORMHR: IPW2,MEM(IPW2)', IPW2, MEM(IPW2)
335  IF( debug ) WRITE(*,*) '% #', iam, ': == pdormhr =='
336  CALL psormhr( 'L', 'N', n, n, ilo, ihi, mem(ipa), 1, 1,
337  $ desca, mem(ipw1), mem(ipq), 1, 1, descq, mem(ipw2),
338  $ -1, info )
339  IF (dpalloc-ipw2.LT.mem(ipw2)) THEN
340  WRITE(*,*) "% Not enough memory for PSORMHR"
341  GO TO 999
342  END IF
343  CALL psormhr( 'L', 'N', n, n, ilo, ihi, mem(ipa), 1, 1,
344  $ desca, mem(ipw1), mem(ipq), 1, 1, descq, mem(ipw2),
345  $ dpalloc-ipw2, info )
346  IF ( info.NE.0 ) THEN
347  WRITE(*,*) "% PSORMHR failed, INFO =", info
348  GO TO 999
349  END IF
350 *
351 * Extract the upper Hessenberg part of A.
352 *
353  CALL pslaset( 'Lower triangular', n-2, n-2, zero, zero,
354  $ mem(ipa), 3, 1, desca )
355 *
356 * Print reduced matrix A in debugging mode.
357 *
358  IF( prn ) THEN
359  CALL pslaprnt( n, n, mem(ipa), 1, 1, desca, 0, 0, 'H', nout,
360  $ mem(ipw1) )
361  CALL pslaprnt( n, n, mem(ipq), 1, 1, descq, 0, 0, 'Q', nout,
362  $ mem(ipw1) )
363  END IF
364 *
365  30 CONTINUE
366  t_hs = mpi_wtime() - t_hs
367 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
368 c $ ' %%% Hessenberg took in seconds:',T_HS
369 *
370 * Compute the real Schur form of the Hessenberg matrix A.
371 *
372  IF( barr ) CALL blacs_barrier(ictxt, 'A')
373  t_qr = mpi_wtime()
374  IF( solver.EQ.1 ) THEN
375  IF( debug ) WRITE(*,*) '% #', iam, ': == pdlaqr1 =='
376 * PRINT*, '% PSLAQR1: IPW1,MEM(IPW1)', IPW1, MEM(IPW1)
377  CALL pslaqr1( .true., .true., n, ilo, ihi, mem(ipa), desca,
378  $ mem(wr1), mem(wi1), ilo, ihi, mem(ipq), descq,
379  $ mem(ipw1), -1, imem, -1, info )
380  IF (dpalloc-ipw1.LT.mem(ipw1)) THEN
381  WRITE(*,*) "% Not enough DP memory for PSLAQR1"
382  GO TO 999
383  END IF
384  IF (intallc.LT.imem(1)) THEN
385  WRITE(*,*) "% Not enough INT memory for PSLAQR1"
386  GO TO 999
387  END IF
388  CALL pslaqr1( .true., .true., n, ilo, ihi, mem(ipa), desca,
389  $ mem(wr1), mem(wi1), ilo, ihi, mem(ipq), descq,
390  $ mem(ipw1), dpalloc-ipw1+1, imem, intallc, info )
391  IF (info.NE.0) THEN
392  WRITE(*,*) "% PSLAQR1: INFO =", info
393  END IF
394  ELSEIF( solver.EQ.2 ) THEN
395  IF( debug ) WRITE(*,*) '% #', iam, ': == pdhseqr =='
396 * PRINT*, '% PSHSEQR: IPW1,MEM(IPW1)', IPW1, MEM(IPW1)
397  IF( barr ) CALL blacs_barrier(ictxt, 'A')
398  CALL pshseqr( 'Schur', 'Vectors', n, ilo, ihi, mem(ipa),
399  $ desca, mem(wr2), mem(wi2), mem(ipq), descq, mem(ipw1),
400  $ -1, imem, -1, info )
401  IF( barr ) CALL blacs_barrier(ictxt, 'A')
402  IF (dpalloc-ipw1.LT.mem(ipw1)) THEN
403  WRITE(*,*) "% Not enough DP memory for PSHSEQR"
404  GO TO 999
405  END IF
406  IF (intallc.LT.imem(1)) THEN
407  WRITE(*,*) "% Not enough INT memory for PSHSEQR"
408  GO TO 999
409  END IF
410  IF( barr ) CALL blacs_barrier(ictxt, 'A')
411  CALL pshseqr( 'Schur', 'Vectors', n, ilo, ihi, mem(ipa),
412  $ desca, mem(wr2), mem(wi2), mem(ipq), descq, mem(ipw1),
413  $ dpalloc-ipw1+1, imem, intallc, info )
414  IF (info.NE.0) THEN
415  WRITE(*,*) "% PSHSEQR: INFO =", info
416  END IF
417  ELSE
418  WRITE(*,*) '% ERROR: Illegal SOLVER number!'
419  GO TO 999
420  END IF
421  t_qr = mpi_wtime() - t_qr
422 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
423 c $ ' %%% QR-algorithm took in seconds:',T_QR
424  t_sch = t_sch + t_qr + t_hs + t_ba
425 * TOTIT = IMEM(1)
426 * SWEEPS = IMEM(2)
427 * TOTNS = IMEM(3)
428  itpereig = float(totit) / float(n)
429  swpspeig = float(sweeps) / float(n)
430  nspeig = float(totns) / float(n)
431 *
432 * Print reduced matrix A in debugging mode.
433 *
434  IF( prn ) THEN
435  CALL pslaprnt( n, n, mem(ipa), 1, 1, desca, 0, 0, 'T',
436  $ nout, mem(ipw1) )
437  CALL pslaprnt( n, n, mem(ipq), 1, 1, descq, 0, 0, 'Z',
438  $ nout, mem(ipw1) )
439  END IF
440 *
441 * Check that returned Schur form is really a quasi-triangular
442 * matrix.
443 *
444  hess = 0
445  DO i = 1, n-1
446  IF( i.GT.1 ) THEN
447  CALL pselget( 'All', '1-Tree', elem1, mem(ipa), i, i-1,
448  $ desca )
449  ELSE
450  elem1 = zero
451  END IF
452  CALL pselget( 'All', '1-Tree', elem2, mem(ipa), i+1, i,
453  $ desca )
454  IF( i.LT.n-1 ) THEN
455  CALL pselget( 'All', '1-Tree', elem3, mem(ipa), i+2, i+1,
456  $ desca )
457  ELSE
458  elem3 = zero
459  END IF
460  IF( elem2.NE.zero .AND. abs(elem1)+abs(elem2)+abs(elem3).GT.
461  $ abs(elem2) ) hess = hess + 1
462  END DO
463 *
464 * Compute residual norms and other results:
465 *
466 * 1) RNORM = || T - Q'*A*Q ||_F / ||A||_F
467 * 2) ORTH = MAX( || I - Q'*Q ||_F, || I - Q*Q' ||_F ) /
468 * (epsilon*N)
469 *
470  stamp = mpi_wtime()
471  IF( compresi ) THEN
472  IF( debug ) WRITE(*,*) '% #', iam, ': Compute residuals 1'
473  IF( debug ) WRITE(*,*) '% #', iam, ': pdgemm 3'
474  CALL psgemm( 'N', 'N', n, n, n, one, mem(ipacpy), 1, 1,
475  $ desca, mem(ipq), 1, 1, descq, zero, mem(ipw1), 1, 1,
476  $ desca )
477  IF( debug ) WRITE(*,*) '% #', iam, ': pdgemm 4'
478  IF( debug ) WRITE(*,*) '% #', iam, ': N=',n
479  IF( debug ) WRITE(*,*) '% #', iam, ': DESCA=',desca(1:dlen_)
480  IF( debug ) WRITE(*,*) '% #', iam, ': DESCQ=',descq(1:dlen_)
481  CALL psgemm( 'T', 'N', n, n, n, -one, mem(ipq), 1, 1,
482  $ descq, mem(ipw1), 1, 1, desca, one, mem(ipa), 1, 1,
483  $ desca )
484  r1 = pslange( 'Frobenius', n, n, mem(ipa), 1, 1, desca,
485  $ dpdum )
486  anorm = pslange( 'Frobenius', n, n, mem(ipacpy), 1, 1,
487  $ desca, dpdum )
488  IF( anorm.GT.zero )THEN
489  rnorm = r1 / (anorm*eps*sqrt(float(n)))
490  ELSE
491  rnorm = r1
492  END IF
493  ELSE
494  rnorm = 0.0d0
495  END IF
496 *
497  IF( comporth ) THEN
498  IF( debug ) WRITE(*,*) '% #', iam, ': Compute residuals 2'
499  CALL pslaset( 'All', n, n, zero, one, mem(ipw1), 1, 1,
500  $ descq )
501  CALL pslacpy( 'All', n, n, mem(ipq), 1, 1, descq, mem(ipw2),
502  $ 1, 1, descq )
503  CALL psgemm( 'T', 'N', n, n, n, -one, mem(ipq), 1, 1, descq,
504  $ mem(ipw2), 1, 1, descq, one, mem(ipw1), 1, 1, descq )
505  o1 = pslange( 'Frobenius', n, n, mem(ipw1), 1, 1, descq,
506  $ dpdum )
507  CALL pslaset( 'All', n, n, zero, one, mem(ipw1), 1, 1,
508  $ descq )
509  CALL psgemm( 'N', 'T', n, n, n, -one, mem(ipq), 1, 1, descq,
510  $ mem(ipw2), 1, 1, descq, one, mem(ipw1), 1, 1, descq )
511  o2 = pslange( 'Frobenius', n, n, mem(ipw1), 1, 1, descq,
512  $ dpdum )
513  orth = max(o1,o2) / (eps*float(n))
514  ELSE
515  orth = 0.0d0
516  END IF
517 *
518  t_res = t_res + mpi_wtime() - stamp
519 c IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
520 c $ ' %%% Residuals took in seconds:',T_RES
521  tottime = mpi_wtime() - tottime
522 c IF( IAM.EQ.0 ) WRITE(*,*)
523 c $ ' %%% Total execution time in seconds:', TOTTIME
524 *
525 *
526 * Print results to screen.
527 *
528  IF( (orth.GT.thresh).OR.(rnorm.GT.thresh) ) THEN
529  passed = 'FAILED'
530  ELSE
531  passed = 'PASSED'
532  END IF
533  IF( debug ) WRITE(*,*) '% #', iam, ': Print results...'
534  IF( iam.EQ.0 ) THEN
535  WRITE( nout, fmt = 9993 ) n, nb, nprow, npcol, t_qr, passed
536  END IF
537  CALL blacs_barrier( ictxt, 'All' )
538  END DO
539  END DO
540  END DO
541  999 CONTINUE
542 *
543 * Deallocate MEM and IMEM.
544 *
545  DEALLOCATE( mem, imem )
546 *
547  CALL blacs_gridexit( ictxt )
548 *
549  777 CONTINUE
550 *
551  CALL blacs_exit( 0 )
552 *
553 * Format specifications.
554 *
555  6666 FORMAT(a2,a3,a6,a4,a5,a6,a3,a3,a3,a9,a9,a9,a8,a8,a9,a9,a9,a9,a9,
556  $ a9,a9,a9,a9,a9,a9,a5,a5,a8,a5,a5)
557  7777 FORMAT(a2,i3,i6,i4,i5,i6,i3,i3,i3,f9.2,f9.2,f9.2,f8.2,f8.2,f9.2,
558  $ f9.2,f9.2,f9.2,f9.2,f9.2,f9.2,f9.2,e9.2,e9.2,e9.2,i5,i5,
559  $ f8.4,i5,i5,a2)
560  9995 FORMAT( ' N NB P Q QR Time CHECK' )
561  9994 FORMAT( '----- --- ---- ---- -------- ------' )
562  9993 FORMAT( i5, 1x, i3, 1x, i4, 1x, i4, 1x, f8.2, 1x, a6 )
563 
564 *
565  END
psgehrd
subroutine psgehrd(N, ILO, IHI, A, IA, JA, DESCA, TAU, WORK, LWORK, INFO)
Definition: psgehrd.f:3
pslaqr1
recursive subroutine pslaqr1(WANTT, WANTZ, N, ILO, IHI, A, DESCA, WR, WI, ILOZ, IHIZ, Z, DESCZ, WORK, LWORK, IWORK, ILWORK, INFO)
Definition: pslaqr1.f:5
pslamch
real function pslamch(ICTXT, CMACH)
Definition: pcblastst.f:7455
max
#define max(A, B)
Definition: pcgemr.c:180
psgebal
subroutine psgebal(JOB, N, A, DESCA, ILO, IHI, SCALE, INFO)
Definition: psgebal.f:2
pslange
real function pslange(NORM, M, N, A, IA, JA, DESCA, WORK)
Definition: pslange.f:3
pslaprnt
subroutine pslaprnt(M, N, A, IA, JA, DESCA, IRPRNT, ICPRNT, CMATNM, NOUT, WORK)
Definition: pslaprnt.f:3
pslacpy
subroutine pslacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pslacpy.f:3
psormhr
subroutine psormhr(SIDE, TRANS, M, N, ILO, IHI, A, IA, JA, DESCA, TAU, C, IC, JC, DESCC, WORK, LWORK, INFO)
Definition: psormhr.f:3
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
pselget
subroutine pselget(SCOPE, TOP, ALPHA, A, IA, JA, DESCA)
Definition: pselget.f:2
pshseqr
subroutine pshseqr(JOB, COMPZ, N, ILO, IHI, H, DESCH, WR, WI, Z, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pshseqr.f:3
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
pshseqrdriver
program pshseqrdriver
Definition: pshseqrdriver.f:8
psmatgen2
subroutine psmatgen2(ICTXT, AFORM, DIAG, M, N, MB, NB, A, LDA, IAROW, IACOL, ISEED, IROFF, IRNUM, ICOFF, ICNUM, MYROW, MYCOL, NPROW, NPCOL)
Definition: psmatgen2.f:4
pslaset
subroutine pslaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: psblastst.f:6863
min
#define min(A, B)
Definition: pcgemr.c:181