R devel - Jan 2006 - bug in qr.coef() and (therefore) in qr.solve (PR#8476)

[I thought I'd submitted this bug report some time ago, but it's never
showed up on the bug tracking system, so I'm submitting again.]

qr.solve() gives incorrect results when dealing with complex matrices or with qr
objects that have been computed with LAPACK=TRUE, whenever the b argument has
more than one column.  This bug flows from qr.coef(), which has a similar
problem.  I believe the problem is the line

  coef[qr$pivot, ] <- .Call("qr_coef_cmplx", qr, y, PACKAGE =
"base")[1:p]

and the similar line in the LAPACK section of the qr.coef() code.  As far as I
can see, this line should read

	coef <- .Call("qr_coef_cmplx", qr, y, PACKAGE =
"base")[qr$pivot,]

With this change, qr.coef() gives correct results for my examples.  In the
examples, the qr.coeffCS() function is my version of qr.coef, with the two lines
in question changed as above (and no other modifications).

Examples:
> A <- matrix(rnorm(9),3,3)
> B <- matrix(rnorm(9),3,3)
> solve(A+1i*B,A+1i*B)
                            [,1]                        [,2] [,3]
[1,]  1.000000e+00+0.000000e+00i -1.853360e-17-1.199306e-17i 0+0i
[2,]  2.338819e-17-1.192988e-19i  1.000000e+00+1.155338e-20i 0+0i
[3,] -6.940188e-18+1.120842e-17i  5.188659e-17-3.226848e-17i
1+0i
> qr.solve(A+1i*B,A+1i*B)
                            [,1]                        [,2]
[1,]  1.000000e-00-2.583088e-16i  1.000000e-00-2.583088e-16i
[2,] -1.045057e-16+1.979352e-16i -1.045057e-16+1.979352e-16i
[3,]  3.966684e-16+7.360601e-16i  3.966684e-16+7.360601e-16i
                            [,3]
[1,]  1.000000e-00-2.583088e-16i
[2,] -1.045057e-16+1.979352e-16i
[3,]  3.966684e-16+7.360601e-16i
## Note:  all columns the same, matrix should be the identity
> qr.coef(qr(A+1i*B),A+1i*B)
                            [,1]                        [,2]
[1,]  1.000000e-00-2.583088e-16i  1.000000e-00-2.583088e-16i
[2,] -1.045057e-16+1.979352e-16i -1.045057e-16+1.979352e-16i
[3,]  3.966684e-16+7.360601e-16i  3.966684e-16+7.360601e-16i
                            [,3]
[1,]  1.000000e-00-2.583088e-16i
[2,] -1.045057e-16+1.979352e-16i
[3,]  3.966684e-16+7.360601e-16i
> qr.coeffCS(qr(A+1i*B),A+1i*B)
                            [,1]                        [,2]
[1,]  1.000000e-00-2.583088e-16i -6.306738e-17-8.893088e-17i
[2,] -1.045057e-16+1.979352e-16i  1.000000e-00-1.921467e-17i
[3,]  3.966684e-16+7.360601e-16i  4.149193e-17-1.149122e-16i
                            [,3]
[1,] -7.350360e-17-6.340421e-17i
[2,] -2.685509e-17+3.225516e-17i
[3,]  1.000000e+00-8.587603e-18i
## Note:  correct results from the function with two modified lines, whether
## LAPACK=TRUE or not and whether complex or not.
> qr.coeffCS(qr(A,LAPACK=TRUE),A)
              [,1]         [,2] [,3]
[1,]  1.000000e+00 0.000000e+00    0
[2,]  1.267908e-15 1.000000e+00    0
[3,] -1.889451e-15 4.074224e-17    1
> qr.coef(qr(A,LAPACK=TRUE),A)
              [,1]          [,2]          [,3]
[1,]  1.000000e+00  1.000000e+00  1.000000e+00
[2,]  1.267908e-15  1.267908e-15  1.267908e-15
[3,] -1.889451e-15 -1.889451e-15 -1.889451e-15
> qr.solve(qr(A,LAPACK=TRUE),A)
              [,1]          [,2]          [,3]
[1,]  1.000000e+00  1.000000e+00  1.000000e+00
[2,]  1.267908e-15  1.267908e-15  1.267908e-15
[3,] -1.889451e-15 -1.889451e-15 -1.889451e-15
> solve(A,A)
             [,1]          [,2] [,3]
[1,] 1.000000e+00  2.628787e-17    0
[2,] 1.899954e-17  1.000000e+00    0
[3,] 2.851722e-16 -6.860465e-17    1
> C <- solve(A,B)
> qr.solve(A,B)-C
              [,1]          [,2]          [,3]
[1,] -8.881784e-16  1.332268e-15  4.440892e-16
[2,] -1.387779e-15  2.220446e-15  8.881784e-16
[3,]  2.664535e-15 -3.552714e-15 -1.110223e-15
## qr.solve() matches solve() with real, non-LAPACK argument
> qr.coef(qr(A),B)-C
              [,1]          [,2]          [,3]
[1,] -8.881784e-16  1.332268e-15  4.440892e-16
[2,] -1.387779e-15  2.220446e-15  8.881784e-16
[3,]  2.664535e-15 -3.552714e-15 -1.110223e-15
> qr.coef(qr(A,LAPACK=TRUE),B)-C
              [,1]      [,2]       [,3]
[1,]  1.110223e-15  3.380377  1.4697337
[2,]  2.164935e-15  1.812489  0.9821958
[3,] -3.552714e-15 -9.280706 -6.3293155
## qr.coef() gives different results with LAPACK=TRUE 
> qr.coeffCS(qr(A,LAPACK=TRUE),B)-C
              [,1]          [,2]          [,3]
[1,]  1.110223e-15 -1.776357e-15 -3.330669e-16
[2,]  2.164935e-15 -3.996803e-15 -6.661338e-16
[3,] -3.552714e-15  7.105427e-15  1.221245e-15
## the modified function gives the same results for LAPACK=TRUE as does
## qr.coef() in the LAPACK=FALSE case.

## lines below just show that qr.coeffCS() works in non-square cases and
## that the problem is there in ths same form in these cases for the original
## qr.coef()
> X <- matrix(rnorm(36),12,3)
> y <- matrix(rnorm(24),12,2)
> b <- qr.coef(qr(X),y)
> qr.coef(qr(X,LAPACK=TRUE),y)-b
              [,1]      [,2]
[1,] -5.551115e-17 0.2509164
[2,]  1.110223e-16 0.1846802
[3,]  0.000000e+00 1.0224349
> qr.coef(qr(X+0i),y)-b
                 [,1]         [,2]
[1,] -5.551115e-17+0i 0.2509164+0i
[2,]  1.110223e-16+0i 0.1846802+0i
[3,]  0.000000e+00+0i 1.0224349+0i
> qr.coeffCS(qr(X+0i),y)-b
                 [,1]            [,2]
[1,] -5.551115e-17+0i 2.775558e-17+0i
[2,]  1.110223e-16+0i 1.110223e-16+0i
[3,]  0.000000e+00+0i 5.551115e-17+0i