Bug: dgelsd() hangs because dlascl divides by zero [solved]

[wdobler]

[I sent this bug report to lapack@cs.utk.edu a week ago, but didn't hear back. Maybe someone on the forum has encountered this before, or even has a fix.]

We are mostly happy users of dgelsd() for fitting data to a model, but recently I found that the routine hangs for certain input data.

What happens is that in dlals0() the variable diflj is zero and then cfrom1 in dlascl() is Inf, which means that that routine's exit criterion is never met.
The matrix has lots of zero entries, but then that is exactly what SVD should be good at.

You find a short test program with a largish (1MB) data file at http://www.capca.ucalgary.ca/~wdobler/tmp/lapack/lapack-bug.tar.gz


Version: 3.1.1 (downloaded two weeks ago from http://www.netlib.org/lapack/lapack.tgz)
OS: Ubuntu 6.06
Compiler: g95, gfortran, g77
CPU: Intel 32-bit, AMD64-bit

Note: This could be related to the bug discussed at http://www.mathkb.com/Uwe/Forum.aspx/num-analysis/3324/SVD-Divide-and-Conquer-LAPACK-DGELSD-problem


Wolfgang Dobler
JPK Instruments AG, Berlin

[wdobler]

Before I start: I am disappointed by the lack of interest from the maintainers' side (if they do exist at all).

Here is a short discussion, plus three patches that fix the problem.

What happened was that the root-finding routine dlasd4() does not always converge within maxit=20 iterations (it needed 23 iterations for the matrix from my bug report), which made dlasd8() return with INFO=1, but dlasd6() ignored this error status and so the call to dgelsd() tried to revert some linear transformations that had actually never been applied (and in the process, division by zero occurred).

The fix consists of three patches:

1. Increase maxit in dlasd4.f from 20 to 200. If you think 200 is a lot, please keep in mind that this will never increase the run time for cases that have previously been successful. dgelsd() will simply try harder in those cases that previously failed.

Code: Select all

--- SRC/dlasd4.f.orig   2006-11-12 22:55:56.000000000 +0100
+++ SRC/dlasd4.f   2007-11-08 16:12:54.000000000 +0100
@@ -94,7 +94,7 @@
 *
 *     .. Parameters ..
       INTEGER            MAXIT
-      PARAMETER          ( MAXIT = 20 )
+      PARAMETER          ( MAXIT = 200 )
       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
      $                   THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0,



2. Fix error handling in dlasd6().

Code: Select all

--- SRC/dlasd6.f.orig   2007-11-12 12:41:27.000000000 +0100
+++ SRC/dlasd6.f        2007-11-12 12:41:27.000000000 +0100
@@ -281,6 +281,13 @@
       CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM,
      $             WORK( ISIGMA ), WORK( IW ), INFO )
 *
+*     Do _not_ try to go on if there was a problem in or below DLASD8,
+*     which could lead to division by zero later.
+*
+      IF (INFO .NE. 0) THEN
+         RETURN
+      ENDIF
+*
 *     Save the poles if ICOMPQ = 1.
 *
       IF( ICOMPQ.EQ.1 ) THEN


3. Print a helpful warning in dlasd8():

Code: Select all

--- SRC/dlasd8.f.orig   2007-11-12 12:41:27.000000000 +0100
+++ SRC/dlasd8.f   2007-11-12 12:41:27.000000000 +0100
@@ -191,6 +191,10 @@
 *        If the root finder fails, the computation is terminated.
 *
          IF( INFO.NE.0 ) THEN
+            WRITE(0,*) 'LAPACK problem: Root finder DLASD4() failed.'
+            WRITE(0,*) 'You can try to get your code working by'
+            WRITE(0,*) 'increasing the parameter MAXIT in dlasd4.f'
+            WRITE(0,*) 'and recompiling the lapack library.'
             RETURN
          END IF
          WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J )


[edrumwri]

I have encountered what appears to be the same bug, but which evokes different behavior. Rather than the routine hanging, LAPACK exits with the message:

** On entry to DLASCL parameter number 4 had an illegal value

Debugging indicates that, just like in Wolfgang Dobler's analysis, diflj in DLALS0 is zero. The result of the multiple divisions (and elements of WORK) then ends up being Inf. My DNRM2 routine returns NaN, instead of Inf, when executed on the WORK array, which is where my results diverge from Wolfgang's. cfrom1 in DLASCL is then NaN, hence the error message that I get above.

I have included a small C++ program that reproduces the behavior on two different machines that I have tested at the end of this message.

Thanks and regards,
Evan

Version: 3.1.2 (Installed from Ubuntu's repository, manually compiled LAPACK yields same behavior)
OS: Ubuntu 10.04 LTS
CPU: Intel 32-bit


C++ test case begins below:


#include <cmath>
#include <algorithm>
#include <assert.h>
#include <string.h>

typedef int INTEGER;
typedef double DOUBLE;
typedef long int __CLPK_ftnlen;

extern "C" {

INTEGER ilaenv_(INTEGER *ispec, char *name__, char *opts, INTEGER *n1,
INTEGER *n2, INTEGER *n3, INTEGER *n4, __CLPK_ftnlen name_len, __CLPK_ftnlen
opts_len);

int dgelsd_(INTEGER *m, INTEGER *n, INTEGER *nrhs,
DOUBLE *a, INTEGER *lda, DOUBLE *b, INTEGER *ldb, DOUBLE *
s, DOUBLE *rcond, INTEGER *rank, DOUBLE *work, INTEGER *lwork,
INTEGER *iwork, INTEGER *info);
}

double log2(double x) { return std::log(x)/std::log(2.0); }

double AA[] = { 0.99999677836890377147, 1.0510329684354502433e-07, 6.1639582327188691124e-13, -4.1211897733983704484e-05, 1.5767387395726473187e-12, 0, -7.5008887989724826184e-12, 0, 0, -8.6068950855278103518e-09, 0, 1.5543122344752191566e-12, -0.00039788540988067166992, -0.00025297374120902826888, -0.00017608123765777694558, -1.8659795930631162264e-11, -0.0007339004748290745006, 0, -4.1009695639360188579e-11, -7.3262507171989454946e-11, -2.6813939957293086991e-11, 0, -1.1544848410593999688e-10, 4.9671378121729503619e-13, 3.5565994593866889772e-13, 0, 0, 8.6262189508135023175e-07, 8.6608707261515149867e-08, -0.00098009818567940110157, -0.00043105336555104845075, -4.556582888781690599e-11, -3.4184877151233195036e-12, -9.9979469148081534513e-11, 8.0201401075896683324e-12, 1.7050366674098427211e-10, 1.355160983429470889e-10, 0, 0, 1.1787348874747749505e-11, 0, 0, 0, 1.0510329684354502433e-07, 0.99999959990478726191, 0, 2.6262325647508077964e-13, 0, 5.5094817597023393319e-13, -7.7627596559715072999e-05, 1.2756462552943048649e-12, 1.3129497489217101247e-12, 0, -6.3896491431236768221e-09, 0, 0, 0, 0, -0.0001775375183865057771, 0, 1.1443068714811488462e-12, -0.00039018518159927850775, -0.00069705334156794362954, -0.00025512036098246371196, 4.2454928461665986106e-13, -0.0010984297064065562388, 0, 0, 3.2757130341565243725e-13, 0, 5.9248841610681779457e-08, 6.8578930645379898579e-08, 0, 0, -0.00043353458581502080094, 0, -0.00095125050066724448072, 1.2166379015354777948e-12, 1.0423995000508057274e-11, 8.2849838101139994251e-12, 1.1367740082590671591e-11, 1.33018596137901568e-11, -1.0355050150678835053e-12, 0, 0, 0, 6.1639582327188691124e-13, 0, 1.0000861725715499695, 0, 8.9046232917278089758e-05, 0, 5.4618113365378917479e-09, 0, 0, -7.1346810174344188482e-05, 0, 0.00039756151784064996946, -2.9744323615243928316e-08, -6.0453215263545079239e-05, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00057798767948596951882, 2.9520153543849403377e-09, 0, 0, -1.7078269489656516811e-05, 0, -7.3257179489250034976e-08, 0, 0, -2.8374846638978823421e-08, 0, -2.643559882109514092e-08, 0, 0, 0, 0, 0, 0, 0, 0, -4.1211897733983704484e-05, 2.6262325647508077964e-13, 0, 0.99996450321456697985, 0, 3.535225112649698076e-09, 0, 0.00010004130162355462375, 2.6131849684141172929e-08, 0, -4.0891095231143648192e-05, 0, 1.6397665669742877981e-08, 1.0425561691729257063e-08, 9.2829661092141968481e-05, 0, -0.00077194398012847553758, 7.3422053992189262317e-09, 0, 0, 1.822845363541603092e-08, 2.7238840072385528401e-09, 0, 0, 0, 2.1017989015881255455e-09, 0.00046979257318469080928, 0, -9.7587177801394808796e-06, 4.0391836320097240787e-08, 0.00022723399845847591649, 0, 0, 0, 0, 0, 0, 0, 0, -1.528752163748592352e-08, 0, 0, 0, 1.5767387395726473187e-12, 0, 8.9046232917278089758e-05, 0, 1.0003162838676038504, 0, 0.00015741992996742837363, 1.9063811640407379855e-08, 1.9712075971511922035e-08, 0.00086733141525546786355, 1.1821213014018283616e-07, 3.3002485133604153589e-08, 0, 0.00084219715887445811831, 0, 0, 0, 0, 0, 1.3452392172430549522e-07, 0, 0, 2.1084096823997811043e-07, 1.0546256756782668162e-08, 7.5513124087045468968e-09, 0, 0, -4.3686543308385505924e-05, 0, 0, 0, 0, -7.2583405952286739193e-08, 0, -6.7622772548237009005e-08, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.5094817597023393319e-13, 0, 3.535225112649698076e-09, 0, 1.000102281124242154, 0, 3.5669418196881963468e-08, 0.00010782931953989827178, 0, -8.5900504189884152595e-05, 0, 0, 0, 0, 0, 0, 0.00039614433692419792621, -3.5367329176416717473e-08, 0, -7.2654553686324341299e-05, 0.00032137844440588425599, 0, 0, 0, 0.00057270806914522243147, 0, 0, -2.0474349611721009978e-05, 0, 0, 0, 0, -8.6207677607585964097e-08, 0, 0, 0, 0, 0, -3.2074097155998515518e-08, 0, 0, 0, -7.5008887989724826184e-12, -7.7627596559715072999e-05, 5.4618113365378917479e-09, 0, 0.00015741992996742837363, 0, 1.0002640628277796875, -3.1608380357539545003e-09, -3.2683222794815947054e-09, 0.00097072711267381928124, -1.9599931455616115272e-08, 1.3773171658737481948e-08, 0, 3.3282297695613749511e-08, 0, 1.3781816354807574498e-08, 0, 0, 3.0289149965945227905e-08, 0.00085489439313518600372, 1.9804388629030000857e-08, 0, 0.0013398867741475917192, 4.4013474531645613297e-09, 3.1514451603875670571e-09, 0, 0, -1.8232028824394230782e-05, 0, 0, 0, 3.3654261433202492526e-08, -3.0291770591883704355e-08, 7.3843319625677139584e-08, -2.8221512737847120889e-08, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.2756462552943048649e-12, 0, 0.00010004130162355462375, 1.9063811640407379855e-08, 3.5669418196881963468e-08, -3.1608380357539545003e-09, 1.0003433966481334316, 0.0003487756047990675512, 1.262095110643457474e-07, 0.00086772496271736443418, 0, 0, 0, 0, 0, 0.0008865043886427947939, 3.5662494124455434985e-08, 0, 0, 2.3596528675451367008e-07, 1.3230424811272456509e-08, 0, 0, 0, 1.0208838663317010287e-08, 0, 0, -4.7399956062232906362e-05, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -7.425441223007922531e-08, 1, 0, 0, 0, 1.3129497489217101247e-12, 0, 2.6131849684141172929e-08, 1.9712075971511922035e-08, 0.00010782931953989827178, -3.2683222794815947054e-09, 0.0003487756047990675512, 1.0003540262202386124, 1.3050126190794486547e-07, 0.00086354355913126346067, 0, 0, 0, 0, 0, 1.5692310340043036376e-07, 3.6709532635903485698e-08, 0, 0, 0.00085932613869044160992, 1.3618865535303115166e-08, 0, 0, 0, 1.0508566627986937192e-08, 0, 0, -4.879160239956714662e-05, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -7.6434496121269290825e-08, 0, 0, 0, -8.6068950855278103518e-09, 0, -7.1346810174344188482e-05, 0, 0.00086733141525546786355, 0, 0.00097072711267381928124, 1.262095110643457474e-07, 1.3050126190794486547e-07, 0.99756598855180866892, 7.8260819547093518622e-07, -0.00018015246809333085309, 0, -0.00043444719071117354758, 0, 0, 0, 0, 0, 8.9479256215607705371e-07, 0, 0, 1.4024191962724863458e-06, -5.7569427725223665249e-05, -4.1220761548776962968e-05, 0, 0, 0.2384741201625984186, 0, 0, 0, 0, 0.00039621500281961097301, 0, 0.00036913612202854606537, 0, 0, 0, 0, 0, 0, 1, 0, 0, -6.3896491431236768221e-09, 0, -4.0891095231143648192e-05, 1.1821213014018283616e-07, -8.5900504189884152595e-05, -1.9599931455616115272e-08, 0.00086772496271736443418, 0.00086354355913126346067, 7.8260819547093518622e-07, 0.9975200472904774962, 0, 0, 0, 0, 0, 9.4517432058793815486e-07, -0.00017864333526634856852, 0, 0, -0.00044259932687129577289, -6.6274871518778866175e-05, 0, 0, 0, -5.113890757629802053e-05, 0, 0, 0.23743954106986669972, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00037196096844849257934, 0, 0, 1, 1.5543122344752191566e-12, 0, 0.00039756151784064996946, 0, 3.3002485133604153589e-08, 0, 1.3773171658737481948e-08, 0, 0, -0.00018015246809333085309, 0, 0.99906865504837583103, 1.2214900479978751946e-07, 0.0002487735854316275308, 0, 0, 0, 0, 0, 0, 0, 0, 0, -5.0729272776217637642e-05, 7.4441631747745873326e-09, 0, 0, -4.3066653842793733986e-05, 0, 3.0084031121813126219e-07, 0, 0, -7.1553484970454661607e-08, 0, -6.6663240372299981118e-08, 0.00056680920916829524714, 0, 0, 0, 0, 0, 0, 0, -0.00039788540988067166992, 0, -2.9744323615243928316e-08, 1.6397665669742877981e-08, 0, 0, 0, 0, 0, 0, 0, 1.2214900479978751946e-07, 0.99883170184660330371, 0.00049111016347941349736, 3.3993686654465449237e-07, 0, 0.00049203832824973314786, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.99797806481424733427, 8.3215432356009699788e-07, 0, 0, 0, 0, 0.00023062737424606183367, 0.00027142527592954879267, 0, 0, 0, 0, 0, 0, -0.00025297374120902826888, 0, -6.0453215263545079239e-05, 1.0425561691729257063e-08, 0.00084219715887445811831, 0, 3.3282297695613749511e-08, 0, 0, -0.00043444719071117354758, 0, 0.0002487735854316275308, 0.00049111016347941349736, 0.99969469083100137574, 4.454407359633094643e-08, 0, 1.8565814968729554835e-07, 0, 0, 0, 0, 0, 0, 2.5122978530855988311e-08, 1.798851140266677362e-08, 0, 0, -0.00010406878171348443729, 0, 0.0012095533584081108458, 1.0904553549018203285e-07, 0, -1.7290602688513345697e-07, 0, -1.6108895373179166199e-07, 0, 0, 0, 0, 0, 0, 0, 0, -0.00017608123765777694558, 0, 0, 9.2829661092141968481e-05, 0, 0, 0, 0, 0, 0, 0, 0, 3.3993686654465449237e-07, 4.454407359633094643e-08, 0.99906876823608470328, 0, 0.00054955030289732809123, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.725689743087288619e-05, 0, 0, 8.3723002858304695906e-07, 0.99839500797929114917, 0, 0, 0, 0, 0, 0.00051178239242044565316, 0, 0, 0, 0, 0, 0, -1.8659795930631162264e-11, -0.0001775375183865057771, 0, 0, 0, 0, 1.3781816354807574498e-08, 0, 0, 0, 0, 0, 0, 0, 0, 0.9990605142780089043, 0, 0, 3.38482343487100934e-07, 0.00055389314693937707901, 4.5293454431494239998e-08, 0, 9.4647545538006916388e-05, 0, -3.1015075221907384417e-05, 0, 0, 0, 0, 0, 0, 0.99838386490855213218, 0, 8.2506132198512105447e-07, 0, 0, 0, 0.00052678582464665835161, 0, 0, 0, 0, 0, -0.0007339004748290745006, 0, 0, -0.00077194398012847553758, 0, 0, 0, 0.0008865043886427947939, 1.5692310340043036376e-07, 0, 9.4517432058793815486e-07, 0, 0.00049203832824973314786, 1.8565814968729554835e-07, 0.00054955030289732809123, 0, 1.0012415562387682044, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3.6272228504108383618e-07, 0, 0, 0.0012117935209980390532, 0.001345271240701995108, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.1443068714811488462e-12, 0, 7.3422053992189262317e-09, 0, 0.00039614433692419792621, 0, 3.5662494124455434985e-08, 3.6709532635903485698e-08, 0, -0.00017864333526634856852, 0, 0, 0, 0, 0, 0, 0.99905516215196821239, 1.240802225410675419e-07, 0, 0.00025543455272236670339, 1.3916702945593328877e-07, 0, 0, 0, -4.380752228155815331e-05, 0, 0, -4.2522576217007834742e-05, 0, 0, 0, 0, 3.024448855803996139e-07, 0, 0, 0, 0, 0.00056683797440360006448, -6.6613751681376953684e-08, 0, 0, 0, -4.1009695639360188579e-11, -0.00039018518159927850775, 0, 0, 0, -3.5367329176416717473e-08, 3.0289149965945227905e-08, 0, 0, 0, 0, 0, 0, 0, 0, 3.38482343487100934e-07, 0, 1.240802225410675419e-07, 0.99883693686762464736, 0.00048692066103067643823, 0.00048585657283484184887, 0, 8.7969789624997574151e-07, 0, 0, 0, 0, 0, 0, 0, 0, 8.2643994975262202729e-07, 0, 0.9980018911044876706, 0, 0, 0, 0.00026164293415437134271, 0.00022972586947078710296, 0, 0, 0, 0, -7.3262507171989454946e-11, -0.00069705334156794362954, 0, 0, 1.3452392172430549522e-07, 0, 0.00085489439313518600372, 0, 0, 8.9479256215607705371e-07, 0, 0, 0, 0, 0, 0.00055389314693937707901, 0, 0, 0.00048692066103067643823, 1.0011731740980491789, 1.7783257327819157467e-07, 0, 0.00018762582209735922945, 0, -4.0678888263423473859e-07, 0, 0, 0, 0, 0, 0, 0.001352341221423469797, 0, 0.0011868325306731430935, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.6813939957293086991e-11, -0.00025512036098246371196, 0, 1.822845363541603092e-08, 0, -7.2654553686324341299e-05, 1.9804388629030000857e-08, 2.3596528675451367008e-07, 0.00085932613869044160992, 0, -0.00044259932687129577289, 0, 0, 0, 0, 4.5293454431494239998e-08, 0, 0.00025543455272236670339, 0.00048585657283484184887, 1.7783257327819157467e-07, 0.99968123772852046649, 2.9467189355258938122e-08, 2.8023189846404505943e-07, 0, 0, 2.2737424221208613062e-08, 0, 0, -0.00010557056922916308039, 0, 0, 1.1060354554448537101e-07, 0, 0.0011842728689849257862, 0, 0, 0, 0, 0, -1.6538159985923783779e-07, 0, 0, 0, 0, 4.2454928461665986106e-13, 0, 2.7238840072385528401e-09, 0, 0.00032137844440588425599, 0, 1.3230424811272456509e-08, 1.3618865535303115166e-08, 0, -6.6274871518778866175e-05, 0, 0, 0, 0, 0, 0, 1.3916702945593328877e-07, 0, 0, 2.9467189355258938122e-08, 1.000105538951036932, 0, 0, 0, -0.00030733160401641423931, 0, 0, -1.57754459254033641e-05, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.4713028512479695564e-08, 0, 0, 0, -1.1544848410593999688e-10, -0.0010984297064065562388, 0, 0, 2.1084096823997811043e-07, 0, 0.0013398867741475917192, 0, 0, 1.4024191962724863458e-06, 0, 0, 0, 0, 0, 9.4647545538006916388e-05, 0, 0, 8.7969789624997574151e-07, 0.00018762582209735922945, 2.8023189846404505943e-07, 0, 0.99766782050310331975, 0, 0.00054670429129971553905, 0, 0, 0, 0, 0, 0, 0.00023107195367649824291, 0, 2.1444208598286884637e-06, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.9671378121729503619e-13, 0, 0.00057798767948596951882, 0, 1.0546256756782668162e-08, 0, 4.4013474531645613297e-09, 0, 0, -5.7569427725223665249e-05, 0, -5.0729272776217637642e-05, 0, 2.5122978530855988311e-08, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.99906190869943545252, 2.3788528369372841098e-09, 0, 0, -1.3762357196689389838e-05, 0, 0, 0, 0, -2.2865593918286464259e-08, 0, -2.1302869968042870141e-08, -2.8916883942908810923e-08, 0, 0, 0, 0, 0, 0, 0, 3.5565994593866889772e-13, 0, 2.9520153543849403377e-09, 0, 7.5513124087045468968e-09, 0, 3.1514451603875670571e-09, 0, 0, -4.1220761548776962968e-05, 0, 7.4441631747745873326e-09, 0, 1.798851140266677362e-08, 0, -3.1015075221907384417e-05, 0, 0, 0, -4.0678888263423473859e-07, 0, 0, 0.00054670429129971553905, 2.3788528369372841098e-09, 0.99905054440775631353, 0, 0, -9.8540990725037858056e-06, 0, 0, 0, -3.0919819925367963975e-05, -8.3149839197327324314e-07, 0, -1.5253243890178680431e-08, 0, 0, -1.6380814649163966124e-08, 0, 0, 0, 0, 0, 0, 3.2757130341565243725e-13, 0, 2.1017989015881255455e-09, 0, 0.00057270806914522243147, 0, 1.0208838663317010287e-08, 1.0508566627986937192e-08, 0, -5.113890757629802053e-05, 0, 0, 0, 0, 0, 0, -4.380752228155815331e-05, 0, 0, 2.2737424221208613062e-08, -0.00030733160401641423931, 0, 0, 0, 0.99904574542793889158, 0, 0, -1.2172623690798189955e-05, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.4989163849298279274e-08, -1.9069026468976346678e-08, 0, 0, 0, 0, 0, 0, 0.00046979257318469080928, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.725689743087288619e-05, 0, -3.6272228504108383618e-07, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.99915808409241513299, 0, -4.2593199000573456203e-06, 0, -1.7182082148514687958e-05, 0, 0, 0, 0, 0, -8.8623504090890037332e-09, 0, 0, -3.1322620444362847536e-05, 0, 0, 0, 8.6262189508135023175e-07, 5.9248841610681779457e-08, -1.7078269489656516811e-05, 0, -4.3686543308385505924e-05, 0, -1.8232028824394230782e-05, 0, 0, 0.2384741201625984186, 0, -4.3066653842793733986e-05, 0, -0.00010406878171348443729, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.3762357196689389838e-05, -9.8540990725037858056e-06, 0, 0, 0.40419017280355429289, 0.064573169076383774012, 0, 0, 0, 4.9654039885993483949e-05, 0, 0.13420123640783032659, 6.4290449088399981292e-05, 5.1097969660240938339e-05, 5.5985923806400883507e-06, 6.5511590177091960641e-06, 4.7530961764796764157e-06, 0, 0, 0, 8.6608707261515149867e-08, 6.8578930645379898579e-08, 0, -9.7587177801394808796e-06, 0, -2.0474349611721009978e-05, 0, -4.7399956062232906362e-05, -4.879160239956714662e-05, 0, 0.23743954106986669972, 0, 0, 0, 0, 0, 0, -4.2522576217007834742e-05, 0, 0, -0.00010557056922916308039, -1.57754459254033641e-05, 0, 0, 0, -1.2172623690798189955e-05, -4.2593199000573456203e-06, 0.064573169076383774012, 0.40218915498390639396, 0, 0, 0, 0, 0, 3.0804014008456448437e-06, 6.4395109748982548581e-06, 5.1181153823964820049e-06, 5.6820939268453773963e-05, 6.6488678471143813908e-05, 0.13596946955251382194, 0, 0, 0, -0.00098009818567940110157, 0, -7.3257179489250034976e-08, 4.0391836320097240787e-08, 0, 0, 0, 0, 0, 0, 0, 3.0084031121813126219e-07, 0.99797806481424733427, 0.0012095533584081108458, 8.3723002858304695906e-07, 0, 0.0012117935209980390532, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.99712653080650182069, 2.0495117104002424924e-06, 0, 0, 0, 0, 0.00023043005525591464533, 0.00027119305127920600995, 0, 0, 0, 0, 0, 0, -0.00043105336555104845075, 0, 0, 0.00022723399845847591649, 0, 0, 0, 0, 0, 0, 0, 0, 8.3215432356009699788e-07, 1.0904553549018203285e-07, 0.99839500797929114917, 0, 0.001345271240701995108, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.7182082148514687958e-05, 0, 0, 2.0495117104002424924e-06, 0.99772241809396822187, 0, 0, 0, 0, 0, 0.00051143699919103946527, 0, 0, 0, 0, 0, 0, -4.556582888781690599e-11, -0.00043353458581502080094, 0, 0, 0, 0, 3.3654261433202492526e-08, 0, 0, 0, 0, 0, 0, 0, 0, 0.99838386490855213218, 0, 0, 8.2643994975262202729e-07, 0.001352341221423469797, 1.1060354554448537101e-07, 0, 0.00023107195367649824291, 0, -3.0919819925367963975e-05, 0, 0, 0, 0, 0, 0, 0.99770839506246655759, 0, 2.0144732941207443844e-06, 0, 0, 0, 0.00052642877647118835327, 0, 0, 0, 0, 0, -3.4184877151233195036e-12, 0, -2.8374846638978823421e-08, 0, -7.2583405952286739193e-08, 0, -3.0291770591883704355e-08, 0, 0, 0.00039621500281961097301, 0, -7.1553484970454661607e-08, 0, -1.7290602688513345697e-07, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.2865593918286464259e-08, -8.3149839197327324314e-07, 0, 0, 4.9654039885993483949e-05, 0, 0, 0, 0, 0.99988030097878344282, 0, -0.00033562193766123460392, 0, 0, 0, 0, 0, 0, 0, 0, -9.9979469148081534513e-11, -0.00095125050066724448072, 0, 0, 0, -8.6207677607585964097e-08, 7.3843319625677139584e-08, 0, 0, 0, 0, 0, 0, 0, 0, 8.2506132198512105447e-07, 0, 3.024448855803996139e-07, 0.9980018911044876706, 0.0011868325306731430935, 0.0011842728689849257862, 0, 2.1444208598286884637e-06, 0, 0, 0, 0, 0, 0, 0, 0, 2.0144732941207443844e-06, 0, 0.9971688362504318448, 0, 0, 0, 0.00026142396049516625212, 0.00022953360777483444721, 0, 0, 0, 0, 8.0201401075896683324e-12, 1.2166379015354777948e-12, -2.643559882109514092e-08, 0, -6.7622772548237009005e-08, 0, -2.8221512737847120889e-08, 0, 0, 0.00036913612202854606537, 0, -6.6663240372299981118e-08, 0, -1.6108895373179166199e-07, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.1302869968042870141e-08, -1.5253243890178680431e-08, 0, 0, 0.13420123640783032659, 3.0804014008456448437e-06, 0, 0, 0, -0.00033562193766123460392, 0, 0.99923697260933685982, 8.3311435528088395586e-10, 6.6215827265736493246e-10, 4.6611459136869370923e-10, 5.4542115268674251638e-10, 8.3282880591895036559e-11, 0, 0, 0, 1.7050366674098427211e-10, 1.0423995000508057274e-11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00056680920916829524714, 0.00023062737424606183367, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.8916883942908810923e-08, 0, 0, 0, 6.4290449088399981292e-05, 6.4395109748982548581e-06, 0.00023043005525591464533, 0, 0, 0, 0, 8.3311435528088395586e-10, 0.99818035840938823533, 7.274721885996981996e-08, 0, 0, 8.764092784829813354e-10, 0, 0, 0, 1.355160983429470889e-10, 8.2849838101139994251e-12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00027142527592954879267, 0, 0.00051178239242044565316, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -8.8623504090890037332e-09, 5.1097969660240938339e-05, 5.1181153823964820049e-06, 0.00027119305127920600995, 0.00051143699919103946527, 0, 0, 0, 6.6215827265736493246e-10, 7.274721885996981996e-08, 0.99825491822400502784, 0, 0, 6.9656902468295811559e-10, 0, 0, 0, 0, 1.1367740082590671591e-11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00052678582464665835161, 0, 0, 0.00026164293415437134271, 0, 0, 0, 0, 0, -1.6380814649163966124e-08, 0, 0, 5.5985923806400883507e-06, 5.6820939268453773963e-05, 0, 0, 0.00052642877647118835327, 0, 0.00026142396049516625212, 4.6611459136869370923e-10, 0, 0, 0.9982238162577203866, 7.1839271975626672884e-08, 5.4929599757613800648e-10, 0, 0, 0, 0, 1.33018596137901568e-11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00056683797440360006448, 0.00022972586947078710296, 0, 0, 0, 0, 0, 0, -2.4989163849298279274e-08, 0, 6.5511590177091960641e-06, 6.6488678471143813908e-05, 0, 0, 0, 0, 0.00022953360777483444721, 5.4542115268674251638e-10, 0, 0, 7.1839271975626672884e-08, 0.99817364672413155802, 6.4275540445635215292e-10, 0, 0, 0, 1.1787348874747749505e-11, -1.0355050150678835053e-12, 0, -1.528752163748592352e-08, 0, -3.2074097155998515518e-08, 0, -7.425441223007922531e-08, -7.6434496121269290825e-08, 0, 0.00037196096844849257934, 0, 0, 0, 0, 0, 0, -6.6613751681376953684e-08, 0, 0, -1.6538159985923783779e-07, -2.4713028512479695564e-08, 0, 0, 0, -1.9069026468976346678e-08, -3.1322620444362847536e-05, 4.7530961764796764157e-06, 0.13596946955251382194, 0, 0, 0, 0, 0, 8.3282880591895036559e-11, 8.764092784829813354e-10, 6.9656902468295811559e-10, 5.4929599757613800648e-10, 6.4275540445635215292e-10, 0.99922886864740467772, 0, 0, 0, -0, -0, -0, -0, -0, -0, -0, -1, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, 0, 0, 0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -1, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, 0, 0, 0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -1, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, -0, 0, 0, 0 };
double xb[] = { 1.1548767587307017809e-11, 8.5059277063119856106e-12, -9.8905760018687825168e-13, -4.894758895534464865e-12, 4.8194195462822224179e-12, -1.0728217454909137448e-12, 8.5032041059102259756e-12, 4.9275373819271468183e-12, 4.9393086636460998542e-12, -3.7249853417498028749e-09, -3.7783238431792148981e-09, -9.4062586998769681982e-12, -6.2507772780446588513e-09, -1.2399920402019842341e-11, -6.2643467447869584545e-09, -6.298316116204532119e-09, -2.3516809816614465565e-11, -1.0748226018220086267e-11, -6.2754999673958048159e-09, -2.2811318408633886735e-11, -1.2244460616462520333e-11, -3.5313456888973366476e-11, -1.7602724560967666727e-10, -1.0076756990932072373e-12, 7.5439588147660057707e-13, -1.0000300410039901391e-12, -6.939428799919463548e-13, 1.0058504852332142526e-08, 1.0306298105584147923e-08, -6.1988536583737862015e-09, -6.3048148957839710295e-09, -6.2855215054590666365e-09, -2.6943686671461207412e-11, -6.179301959703759736e-09, -2.089849231221382072e-09, -4.85936973384906224e-12, -1.1616117994280668915e-11, -1.1804627042652756697e-11, -4.9024424334135392083e-12, -2.1569519352847458808e-09, 0, 0, 0 };

/// Calls LAPACK function for least-squares solution to Ax=b using SVD
static void gelsd_(INTEGER* M, INTEGER* N, INTEGER* NRHS, DOUBLE* A, INTEGER* LDA, DOUBLE* B, INTEGER* LDB, DOUBLE* RCOND, INTEGER* INFO)
{
// create array to hold singular values
INTEGER min_mn = std::min(*M, *N);
double* S = new double[min_mn];

// necessary for computing WORK sizes
INTEGER ISPEC = (INTEGER) 9;
INTEGER TMP = (INTEGER) 0;
const char* NAME = "DGELSD";
const char* OPTS = " ";
INTEGER smlsiz = ilaenv_(&ISPEC, (char*) NAME, (char*) OPTS, M, N, NRHS, &TMP, strlen(NAME), strlen(OPTS));
assert(smlsiz > (INTEGER) 0);
INTEGER NLVL = std::max((INTEGER) 0, (INTEGER) (log2((double) min_mn/(double) (smlsiz+1))+1));
INTEGER LIWORK = std::max((INTEGER) 1, 3*min_mn*NLVL + 11*min_mn);
INTEGER* IWORK = new INTEGER[LIWORK];
INTEGER RANK;

// first do WORK query
double WORK_QUERY;
INTEGER LWORK = -1;
dgelsd_(M, N, NRHS, A, M, B, LDB, S, RCOND, &RANK, &WORK_QUERY, &LWORK, &LIWORK, INFO);
assert(*INFO == 0);

// setup LWORK
LWORK = (INTEGER) WORK_QUERY;

// setup WORK array
double* WORK = new double[std::max((INTEGER) 1, LWORK)];

for (INTEGER i=0; i< LWORK; i++)
WORK[i] = (double) 0.0;
for (INTEGER i=0; i< LIWORK; i++)
IWORK[i] = (INTEGER) -1;
// compute
dgelsd_(M, N, NRHS, A, M, B, LDB, S, RCOND, &RANK, WORK, &LWORK, IWORK, INFO);

delete [] IWORK;
delete [] S;
delete [] WORK;
}

int main()
{
// setup LAPACK parameters
INTEGER M = 43;
INTEGER N = 43;
INTEGER NRHS = 1;
INTEGER LDA = M;
INTEGER LDB = std::max(M,N);
INTEGER INFO;
double stol = -1.0;

gelsd_(&M, &N, &NRHS, AA, &LDA, xb, &LDB, &stol, &INFO);
}