update 2009 to 2010 in some files

Start getting ready for the release

Update CMake build

fix again some COMPLEX instead of CMPLX

fix again some COMPLEX instead of CMPLX

Fix a minor problem on csd - intrinsic function is CMPLX and not COMPLEX (same for double precision)

Correct bug0064 :: got a bug in divide and conquer (non convergence)
    o see http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=529 (forum topic 529)
    o Increase MAXIT to 200 (was 20 before)
    o Add some XERBLA calls when INFO is not zero.

Incorporated CSD source and testing routines into SRC/Makefile, TESTING/Makefile, TESTING/EIG/Makefile, TESTING/EIG/xchkee.f; changes to comments in CSD routines

Removed nonexistent fourth matrix type from CSD test routines and csd.in

Modified TESTING/EIG/alahdg.f to test CS decomposition code

Added CS decomposition test files to TESTING

Added CS decomposition source files to SRC/

Commiting the 3 other precisions (single, complex, dcomplex) for sytri using Level BLAS 3.
Update testing accordingly

 [xSYTRS/xSYSV] Hide the call to syconv inside trs2 to avoid changing the SYTRS interface.
Update the testing accordingly

 [DSYTRI2]      Comit dsytri2 to get some feedback
                Update the testing accordingly

DSYTRI2 is the Level 3 blas Version of DSYTRI
The actual routine that does the work is DSYTRI2X (name can be changed)
DSYTRI2 is just a wrapper to allow to hide the 2D Workspace required by the routine.
The interface had to be changed to integrate the possibility of doing a workspace query.
DSYTRI2x implementation will be documented in a LAWN. This algorithm was inspired by the following paper:
"Families of Algorithms Related to the Inversion of a Symmetric Positive Definite Matrix"
PAOLO BIENTINESI Duke University and BRIAN GUNTER Delft University of Technology and ROBERT A. VAN DE GEIJN The University of Texas at Austin

Add hetrs2 for c and z, a Level BLAS 3 version of hetrs (same than rev 751)

Included bug fix provided by Zlatco on Jacobi SVD
Email from Zlatco on August 24th 2010:
The problem that was reported (with zero matrix) is caused by bad initialization to xLASSQ.
It should be ZERO, ONE and not ZERO, ZERO. In fact, I had it ZERO, ONE throughout
the complete development of the code and decided to change it to ZERO, ZERO a the very
end to make it "more elegant". That was stupid, because xLASSQ does not touch those
variables in case of zero vector, leaving scaling at ZERO, and in the nonzero case the scaling
is between ONE and SQRT(N). So, in case of zero vector, division by a variable that
 is normally bigger than ONE causes division by zero.
I have corrected that and few other things, stress tested the code and it should be OK now.

README:
i)   In xgejsv.f and xgesvj.f input parameters SCALE and
     SUMSQ in xlassq.f are now initially set as SCALE = ZERO, SUMSQ=ONE.
     Setting them both to zero (without carefully reading xlassq.f) caused
     problems with exactly zero columns.
ii)  There was a problem in the branch that computes only SIGMA and U of a
     rank deficient matrix. The computed numerical rank (NR) was incorrectly
     written as N in parameter lists of the corresponding calls.
iii) In xgsvj0.f, xgsvj1.f testing the input parameters is changed to prevent
     unnecessarily negative INFO in some situations.
iv) Minor changes, renaming some variables etc.

Array out-of-bounds reference in xLAQR5.

Bug report and fix from Mathew Cross (NAG) on Sat Aug 14 2010:

In the section of the code "Special case: 2-by-2 reflection (if needed)" the logical test
IF( BMP22 .AND. ( V( 1, M22 ).NE.ZERO ) )
must be split into a nested test
IF (BMP22) THEN, IF (V(1,M22).NE.ZERO) THEN...
If .NOT. BMP22 then M22 can exceed the second extent of V (and recall that logical expressions can be evaluated in any order in Fortran).

See http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=1949.

Patch provided by the kitware team on Aug 5th (Brad King and Bill Hoffman)

Add the install rules and exports the targets.
Now one can write

find_package(LAPACK 3.2.2)

in another project to find either the LAPACK build tree or the
install tree.

Various comment fixes

Adding a hint on the data storage required for band matrices in xLASCL.

Split some lines (continue)

Split some lines

 Add static/dynamic library, package support and rearrange source files

Time has come to have the fortran90 slamch and dlamch in the lapack package.

Jason (Riedy) wrote our ( ... his? :) ) ideas about it three years ago:
http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289
Piotr (Luszczek) has written two subroutines, tested them on few platforms, collected the result on his webpage, and sent emails to the lapackers a few times. See:
http://www.cs.utk.edu/~luszczek/lapack/lamch.html
Theses slamch.f and dlamch.f subroutines were taken from PLASMA-2.1.0.

Change to the LAPACK library:
* move the current LAPACK subroutine slamch.f (resp dlamch.f) as
 slamchf77.f (resp. dlamchf77.f),
* take the new slamch.f subroutines (resp. dlamch.f), remove the PLASMA
 header, have a LAPACK header, and insert the new routines in the
 library.

Minor:
* I would leave these routines compiled with the NOOPT flag.

Problem:
* CLAPACK: no idea how CLAPACK's going to handle this. CLAPACK can rely on
 IEEE arithmetic, can relay on float.h, or can rely on the previous
 xlamch.f

Corrct bug0060: Typo in TESTING/EIG/dchkee.f and TESTING/EIG/zchkee.f (bug report sent by Ake Sandgren to lapackers on July 1st 2010)

Fix bug0060 ::    Bug in dtgsy2, missing IWORK call
bug report and fix sent by Vasile Sima on 24 Feb 2009 to "lapack@cs.utk.edu".

The error is that "IWORK" must appear in the line 640.  (It was included in the 1999 version.)
Note that the statement IE = ( I+1 ) - 1 is strange;  also, in other similar place, IWORK is used like below.

Fix Bug0059
bug report and fix sent by Alexander Kobotov (Intel) on Wed, 23 June 2010 to "lapack@cs.utk.edu".

ZGELSD crashes on large matrices (eq M=N=NHRS>=88), if real workspace allocated exactly as reported on workspace query.
Root cause is wrong estimation of real workspace size in ZLALSA (on path zgelsd->zlalsd->zlalsa). It assumes that RWORK array length should be max( N, (SMLSZ+1)*NRHS*3 ), where SMLSIZ=25. But in fact the function calls zlals0 which requires RWORK to be ( K*(1+NRHS) + 2*NRHS ), where K could be as large as number of columns in matrix A.
The same is actual for CGELSD.

Update date and version for 3.2.2 release

Corrected a bug in SBDSDC and DBDSDC.

The routine was returning with INFO = 0 while a failure of convergence should
have been reported with INFO > 0.

The bug has been around for quite a while and has been difficult to reproduce.

More information:

(*) bug0025 :: serious problem in DGESDD: DGESDD returns negative singular values with INFO = 0

    o reported by Guy Bencteux on Sat Dec 06 2008
    o see forum topic 847
    o Julie : CANNOT REPRODUCE THE PROBLEM

    o on Sat Jun 19 2010, Bruno Pinçon from ESIAL, reported another matrix with the same problems
    o see forum topic 1868
    o we can reproduce the problem with on MacOS with gfortran reference LAPACK
      and vecLib BLAS, we can not reproduce the problem with reference LAPACK and reference BLAS.
    o A patch has been proposed by Julien Langou. (See forum topic 1868.)
    o Under review by other lapackers.
    o Note: there is strong evidence that the bug of Guy Bencteux and this one
      are the same, although we were never able to reproduce Guy's problem. So
      we assume Guy's problem gone! (And send an email to Guy to inform him.)

Update CMAKE build with new routines

OK. So I got upset and did a
  grep "an singular value" *f
in our SRC files to fix them all.
Here we go.

another typo

typo in the comment

Step 2 of xlarfp: add new routines and add same test code plus check R(i,i) is nonnegative and real

 Step 1 of xlarfp : back out the old changes (need to run testing before going on)

 Missing local variable declaration (IINFO)

 Missing local variable declaration

Add SYTRS2 routine - A BLAS 3 version of SYTRS
Add SYCONV routine: convert back and forth the factorization returned by SYTRF to be able to call SYTRS2.
Modify SYSV that now is calling SYTRS2 instead of SYTRS (and also SYCONV to convert and revert the factorization returned by SYTRF).
Modify testing to have TRS but also TRS2 tested in the LIN testing for SY.

The variable N0 in xLASQ3 is INPUT/OUTPUT. (The comments said INPUT.)

The bug report and bug fix are from William Gandler (NIH/CIT). The related
emails are on the mailing list dated as follows: October 16, 2009; April 23,
2010; May 01, 2010, May 03, 2010.

typo in comments spotted by Jim

.. fixed some comments ..

from comments from A. Buttari, CNRS, France, to Julien Langou on 01/20/2010

remove comment to not include xerbla and xerbla_array

Many comment fix for allowing C Wrapper to work on MATGEN

Add descriptions of the arguments work and uplo following Paul Roberts[NAG]'s comment sent to lapack@cs.utk.edu

Fix comment following Paul Roberts'email on lapack mailing-list
I was just looking at zgbsvxx.f from LAPACK 3.2.1 and I noticed a small
mistake in some of the argument descriptions. AB, AFB, B and X are all
down as DOUBLE PRECISION rather than COMPLEX*16.

Fix issue in testing suite found by Michael Chuvelev

Fix comments in ilaprec and remove some space in Makefile

   Following bug report and fix sent by Michael Chuvelev (Intel) on Fri, 2 Apr 2010 to "lapack@cs.utk.edu".
Issue in LAPACK 3.2.1 *la_{po,sy,he}rfsx_extended.f sources. For instanse, in sla_syrfsx_extended.f:

DY is described as
*     DY             (input) REAL array, dimension (N)
*     Workspace to hold the intermediate solution.

whereas N-by-NRHS array is referenced at DY due to the following statement (line 384):

            CALL SSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, DY, N, INFO )

This may lead to an outside memory corruption. I assume one RHS should be referenced here, because all the RHSs are considered individually in the loop by J. That is, substituting NRHS by 1:

            CALL SSYTRS( UPLO, N, 1, AF, LDAF, IPIV, DY, N, INFO )

should fix the issue. This is actual for any subroutine mentioned.

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

"Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by Jim Demmel
and Guillaume Revy. See forum post 1783.

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Fix problem in xTGEX2. The threshold value is too stringent and some matrices
are failing. Relax the threshold by a factor 2. More below.

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

See: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=1783
dgges issue >> by OndraKamenik >> Tue Mar 16, 2010 6:01 am

All,

I have the following problem with dgges. For version 3.1.1 and sooner, I get a
reasonable result, for version 3.2 and 3.2.1 I get info=n+2.

I am separating eigenvalues in the unit circle from one outside the unique
circle.

The two D and E matrices have relatively well separated null spaces, the
minimum angle is acos(0.97). However, if i calculate condition numbers of
E-lambda*D of lambda=[-1:0.01:1], they are quite bad.

The matrices are attached with a small c++ program which calls dgges and sorts
eigenvalues for a better comparison. There is also a Makefile, which links with
different version of lapack.

I use the reference blas.

My question is if it is a bug in Lapack introduced between 3.2. and 3.1.1 or
the matrix is just very bad and in version 3.1.1 I was just lucky to get a
reasonable solution.

Many thanks for any help.

Ondra Kamenik

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Much more conversation on mailing list and forum [ skipped ]

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Date: Wed, 31 Mar 2010 19:53:11 -0600
From: James Demmel

This bug was introduced by changing the routine that computes Householder
transformations in the last release of LAPACK from dlarfg to the new dlarfp,
which makes the diagonal of the R factor in QR nonnegative. The bug was a
failure (INFO = N+3) in the dgges routine for computing selected deflating
subspaces of a matrix pencil A - lambda*B. I will describe the bug, a quick
fix, and implications for the floating point debugging project that Guillaume
and others of us are working on.

The bug occurred in dtgex2, in the swapping of eigenvalues (adjacent 1x1 and/or
2x2 blocks on the diagonal of the generalized Schur form) in order to compute
a selected subspace. The swapping involves QR decompositions of  small (n <= 4)
matrices.

The code in dtgex2 performs two of its own internal correctness tests (a "weak"
one and a "strong" one) to see if the swapping has been performed stably. The
two tests compute a residual in slightly different ways. The weak test passed,
but the strong test  (which can be commented out by setting the internal
parameter WANDS to be .false.)  failed, leading to returning INFO = N+3.
However, it failed by exceeding the threshold THRESH only by a factor like 1.2
or less. THRESH is set to 10*macheps*dnorm, so if we changed the (probably
somewhat arbitrary) factor from 10 to 20, it would work. Or we could set WANDS
= false.

Since Bo's name is on this routine, his comments are particularly welcome.

Guillaume and I spent a while tracking this down, and will continue to find out
why dlarfp led to a (slightly!) larger residual than the old dlarfg.

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Date: Thu, 1 Apr 2010 20:01:21 -0600
From: James Demmel <demmel@cs.berkeley.edu>

Actually, on reexamining the data, this example would be fixed by
changing 10 to 11,
but let's go with 20 to be on the safe side :) . I think it would be
most efficient
if Julie or you appropriately change the one line of code in DTGEX2:

      THRESH = MAX( TEN*EPS*DNORM, SMLNUM )

All versions (S/D/C/Z) have an analogous line of code with the same
constant TEN that I would change.

Sorry, I've lost track of the Mathworks bug report on dlarfp. Can you
remind me?

Thanks,
Jim

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Fix CMakeLists.txt only blas is exported, lapack needs to be exported as well. (From Bill Hoffman)

Fix comments to generate wrapper correctly

Correct bug0050 :: An extra argument after matrix A lines 348-249 in STFSM.f causes a crash of SGEMM call sent from Alexander Kobotov (Intel) on Wed, 3 Feb 2010 to lapack@cs.utk.edu

bug0049 :: output parameter GIVPTR stays uninitialized in case of quick exit (n.eq.0) or if the rank-1 modifier is small enough in *laed8 sent by Alexander Kobotov (Intel) on Wed, 3 Feb 2010 to lapack@cs.utk.edu

Fix bug0048-Hanging could occur in *gebal if a NaN is in input matrix sent by Alexander Kobotov(Intel) to lapack on Feb 3rd 2010

Fix bug0046-Incorrectly documented RWORK workspace in ZGESDD sent by user Zbigniew on Forum (topic 1779)

Take off comment to put xerbla back in BLAS lib

Fix buggy comments reported by Intel folowing Jim's recommendation

Shorten long lines

Add testversion and testieee as executables

Following Jim and Sven report,  sgeqpf and shseqr are computational routines.  slanv2 is an auxiliary routine

Some minor comments modifications

Add make.inc for ifort

Following a complain from John Tellefson (Salina, KS), added the make of the
variant library to the testing of the variants if the variant library is not
present.

Bug 0046 fixed.

author: Vasile Sima
committer: Julien Langou

Corrected LAPACK routines dlagv2 and slagv2 based on Vasile Sima (National
Institute for Research & Development in Informatics, Bucharest, Romania)'s
email to lapack@cs.utk.edu on Monday 25 January 2010.

> Specifically, the variable WI, which is used in the line 268 (close to the
> end), is not initialized (with ZERO) in the cases "A can be deflated" and "B
> is singular" (i.e., the code segments in the lines 138-144, 148-156, and
> 158-167).  The corrected version is included in the attached archive.

Polish some comments, etc.. for the C wrapper

Remove a trailing blank line at the end of dlarfp.f following Andy May (Cardiff
University) 12/31/09's comment.

Still working on there subroutines ...  I have added a variable SAVEALPHA, I am
not sure whether it is needed or not.  What I am pretty sure of, is that the
code is correct with, the previous code without may be correct as well ...
Anyway, I like it better like this.

Darned!!! I messed up with my commit r709!!!!

This is
            " ( 1 - CONJG( TAU ) ) * ( ALPHA ) = ABS( ALPHA ) "
that needs to be true so, indeed, the previous code was correct and r709 is a
mistake. (r710 is good.) So this commit rollbacks r709 and modifies a little
r710.

This is the relevant change. This is the bug fix for the MathWorks/PatQuillen's
bug in xLARFP. The fix seems to work her in Denver ...

After three irrelevant commit ( r706, r707, r708), this is the first relevant
modifications. I believe that in the complex case, when NORMX = ZERO but ALPHI
is not ZERO (i.e. ALPHA has a complex imaginary part), then TAU needs to be
such that:

            " ( 1 - TAU ) * ( ALPHA ) = ABS( ALPHA ) "

Since we have
            ALPHR = REAL( ALPHA )
            ALPHI = AIMAG( ALPHA )
            XNORM = DLAPY2( ALPHR, ALPHI )
The way to do this is to set TAU with

            TAU = CMPLX( ONE - ALPHR / XNORM, ALPHI / XNORM )

as opposed to

            TAU = CMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )

(Note: XNORM is used as temporary variable here)

change SAFMIN with SMLNUM
change RSAFMN with BIGNUM

remove the comments on the constraints of TAU
this comment is true for xLARFG, it is not true anymore for xLARFP

start some modification of the xLARFP routines
modification 1: remove the comments starting with ! and start them with *

Correct a bug found by Ashutosh Mahajan from Mathematics and Computer Science Division of Argonne National Labs.
Email sent on Thu, 10 Dec 2009 16:53:57 to lapack@cs.utk.edu.

I think we have a bug in dstemr.f for the case when N is 2. The following
block:

392:                   IF (CS.NE.ZERO) THEN
393:                      ISUPPZ(2*M-1) = 1
394:                      ISUPPZ(2*M-1) = 2
395:                   ELSE
396:                      ISUPPZ(2*M-1) = 1
397:                      ISUPPZ(2*M-1) = 1
398:                   END IF

should really be (note lines 394, 397):

392:                   IF (CS.NE.ZERO) THEN
393:                      ISUPPZ(2*M-1) = 1
394:                      ISUPPZ(2*M) = 2
395:                   ELSE
396:                      ISUPPZ(2*M-1) = 1
397:                      ISUPPZ(2*M) = 1
398:                   END IF

similarly the block:
415:                IF (SN.NE.ZERO) THEN
416:                   IF (CS.NE.ZERO) THEN
417:                      ISUPPZ(2*M-1) = 1
418:                      ISUPPZ(2*M-1) = 2
419:                   ELSE
420:                      ISUPPZ(2*M-1) = 1
421:                      ISUPPZ(2*M-1) = 1
422:                   END IF
423:                ELSE

should really be (note lines 418, 421):

415:                IF (SN.NE.ZERO) THEN
416:                   IF (CS.NE.ZERO) THEN
417:                      ISUPPZ(2*M-1) = 1
418:                      ISUPPZ(2*M) = 2
419:                   ELSE
420:                      ISUPPZ(2*M-1) = 1
421:                      ISUPPZ(2*M) = 1
422:                   END IF
423:                ELSE

Correct a bug found by Ashutosh Mahajan from Mathematics and Computer Science Division of Argonne National Labs.
Email sent on Fri, 11 Dec 2009 14:51:07 to lapack@cs.utk.edu.

Subject: [Lapack] bug in dsyevr.f when N = 1

when dsyevr is called with N=1, then it returns without setting up ISUPPZ values.

331:       IF( N.EQ.1 ) THEN
...
342:          IF( WANTZ )
343:      $      Z( 1, 1 ) = ONE
344:          RETURN
345:       END IF

It should rather do

IF( WANTZ )
    Z( 1, 1 ) = ONE
    ISUPPZ(1) = 1
    ISUPPZ(2) = 1
END IF

 some more typos

typo in comments

Fix bug0038 - dgelsd waas not returning iwork size

Bug report from Michael Chuvelev from Intel on DSYEVR.

======================================================================================
Date: Wed, 14 Oct 2009 02:42:06 -0600
From: "Chuvelev, Michael" <michael.chuvelev@intel.com>
To: "lapack@cs.utk.edu" <lapack@cs.utk.edu>
Subject: [Lapack] VL,   VU are referenced in lapack-3.2.1 dsyevr even if range.ne.'V'

Hello,

dsyevr contains this code (lines 360-361):

      VLL = VL
      VUU = VU

which means VL, VU is accessed regardless of RANGE.

Whereas this issue is fixed in ssyevr, for instance (lines 364-367):

      IF (VALEIG) THEN
         VLL = VL
         VUU = VU
      END IF

Best regards,

Michael.
======================================================================================

Correct problem reported by Kevin Wadleigh to keep compilers from optimizing NAN6 = NAN5*0.0

Date: Tue, 29 Sep 2009 06:17:06 -0600
From: "Chuvelev, Michael"
Subject: [Lapack] Bug in lapack-3.2.1/TESTING/EIG/zdrvgbx.f

 Hi,

please pay attention to the printing format 9998 misuse in zdrvgbx.f:

                              WRITE( NOUT, FMT = 9998 )'ZGBSVXX', FACT,
     $                             TRANS, N, KL, KU, IMAT, 7,
     $                             RESULT( 7 )
...

9998 FORMAT( ' *** In ZDRVGB, LAFB=', I5, ' is too small for N=', I5,
     $      ', KU=', I5, ', KL=', I5, /
     $      ' ==> Increase LAFB to at least ', I5 )

This would possibly lead to a segmentation fault in case of errors in zgbsvxx.

The same issue concerns other precisions too.

Best regards,

Michael.

 Fix bug0038 see out of range element access in DLASET when using DGELSS (forum topic 1582)

Commit Julien's proposed change for bug0019

remove extra definition

Fix whitespace comments detected from parser

Fix bug0023: SROTMG and DROTMG uses deprecated Fortran ASSIGN statement and assigned GOTO statement, actually fixed ROTM also

 Fix ifort flags just for Unix

 Fix variables used

 Use try_compile to look for the timing routine

Add conf file for submitting machine to LAPACK Dashboard

Better solution as we set NONE as default...WINDOWS build run great...

Need the NONE for Windows...

Add CMAKE support to LAPACK

Fixed comments on factorization in refinement routines. Found when coding in C++

In refinement routines, err_bnds initialization if-else statements should be if statements, found when coding in C++

thru -> through

fixed the variants_testing Makefile problem

Bug report sent by Alexander V. Kobotov (from Intel) on Mon, 6 Apr 2009 to "lapack@cs.utk.edu".
"c(he/sy)equb: no slamch in external lists and lsame, intrinsics not described"

[[ I did not really want to commit the previous commit, this is the follow-up commit ... ]]

Bug report sent by Alexander V. Kobotov (from Intel) on Mon, 6 Apr 2009 to "lapack@cs.utk.edu".

"(d/s)tgsen.f: iwork(1) always referenced: line 455: iwork( 1 ) = LIWMIN,
documentation says that if IJOB=0 it shouldn't, so NULL pointer causes a
sigfault."

(There is indeed the same problem for the WORK array and the IWORK array in the
complex routines.)

I have changed the header of the routines ctgsen.f, dtgsen.f, stgsen.f, and
ztgsen.f. A Fortran array needs to be of size at least 1. So IWORK is of size
at least 1. It was indeed written in the header of the routine:
"IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))"
So since IWORK is of size at least 1, there is no reason not to reference it.
When there is a workspace query with IJOB.EQ.0, WORK(1) and IWORK(1) are both
set to 1.

Bug report sent by Alexander V. Kobotov (from Intel) on Mon, 6 Apr 2009 to "lapack@cs.utk.edu".

"(d/s)tgsen.f: iwork(1) always referenced: line 455: iwork( 1 ) = LIWMIN,
documentation says that if IJOB=0 it shouldn't, so NULL pointer causes a
sigfault."

There is indeed the same problem for array work.

I have changed the header of the routines ctgsen.f, dtgsen.f, stgsen.f, and
ztgsen.f. No matter what an array needs to be of size at least 1. So
IWORK is of size at least 1. It was written in the header of the routine:
"IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))"
So since IWORK is of size at least 1, there is no reason not to reference it.
When there is a workspace query with IJOB.EQ.0, IWORK(1) is set to 1.