QRsolve(a, b)  
          or QRsolve(a, b, which=which)  
  
  
     returns the solution x (in a least squares or least norm sense  
     described below) of the matrix equation:  
        A(,+)*x(+) = B  
     If A is an m-by-n matrix (i.e.- m equations in n unknowns), then B  
     must have length m, and the returned x will have length n.  
  
     If nm, the system is underdetermined -- many solutions are possible  
             -- the returned x has minimum L2 norm among all solutions  
  
     B may have additional dimensions, in which case the returned x  
     will have the same additional dimensions also have those dimensions.  
     The WHICH dimension of B and the returned x is the one of length m  
     or n which participates in the matrix solve.  By default, WHICH=1,  
     so that the equations being solved are:  
        A(,+)*x(+,..) = B  
     Non-positive WHICH counts from the final dimension (as for the  
     sort and transpose functions), so that WHICH=0 solves:  
        A(,+)*x(..,+) = B  
  
     QRsolve works by QR factorization if nm.  
     QRsolve is slower than LUsolve.  Its main attraction is that it can  
     handle overdetermined or underdetermined systems of equations  
     (nonsquare matrices).  QRsolve may fail for singular systems; try  
     SVsolve in this case.  
  

unknown type function, documented at startup/matrix.i   line 337