------------------------------------------------------------------------------- |-- Module : Numeric.LinearAlgebra.LAPACK-- Copyright : (c) Alberto Ruiz 2006-7-- License : GPL-style-- -- Maintainer : Alberto Ruiz (aruiz at um dot es)-- Stability : provisional-- Portability : portable (uses FFI)---- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>).-------------------------------------------------------------------------------moduleNumeric.LinearAlgebra.LAPACK(-- * Matrix productmultiplyR,multiplyC,multiplyF,multiplyQ,-- * Linear systemslinearSolveR,linearSolveC,lusR,lusC,cholSolveR,cholSolveC,linearSolveLSR,linearSolveLSC,linearSolveSVDR,linearSolveSVDC,-- * SVDsvR,svRd,svC,svCd,svdR,svdRd,svdC,svdCd,thinSVDR,thinSVDRd,thinSVDC,thinSVDCd,rightSVR,rightSVC,leftSVR,leftSVC,-- * EigensystemseigR,eigC,eigS,eigS',eigH,eigH',eigOnlyR,eigOnlyC,eigOnlyS,eigOnlyH,-- * LUluR,luC,-- * CholeskycholS,cholH,mbCholS,mbCholH,-- * QRqrR,qrC,-- * HessenberghessR,hessC,-- * SchurschurR,schurC)whereimportData.Packed.InternalimportData.Packed.MatriximportNumeric.ConversionimportNumeric.GSL.Vector(vectorMapValR,FunCodeSV(Scale))importForeign.Ptr(nullPtr)importForeign.C.TypesimportControl.Monad(when)importSystem.IO.Unsafe(unsafePerformIO)-----------------------------------------------------------------------------------foreignimportccallunsafe"multiplyR"dgemmc::CInt->CInt->TMMMforeignimportccallunsafe"multiplyC"zgemmc::CInt->CInt->TCMCMCMforeignimportccallunsafe"multiplyF"sgemmc::CInt->CInt->TFMFMFMforeignimportccallunsafe"multiplyQ"cgemmc::CInt->CInt->TQMQMQMisTMatrix{order=ColumnMajor}=0isTMatrix{order=RowMajor}=1ttx@Matrix{order=ColumnMajor}=xttx@Matrix{order=RowMajor}=transxmultiplyAuxfstab=unsafePerformIO$dowhen(colsa/=rowsb)$error$"inconsistent dimensions in matrix product "++show(rowsa,colsa)++" x "++show(rowsb,colsb)s<-createMatrixColumnMajor(rowsa)(colsb)app3(f(isTa)(isTb))mat(tta)mat(ttb)matsstreturns-- | Matrix product based on BLAS's /dgemm/.multiplyR::MatrixDouble->MatrixDouble->MatrixDoublemultiplyRab={-# SCC "multiplyR" #-}multiplyAuxdgemmc"dgemmc"ab-- | Matrix product based on BLAS's /zgemm/.multiplyC::Matrix(ComplexDouble)->Matrix(ComplexDouble)->Matrix(ComplexDouble)multiplyCab=multiplyAuxzgemmc"zgemmc"ab-- | Matrix product based on BLAS's /sgemm/.multiplyF::MatrixFloat->MatrixFloat->MatrixFloatmultiplyFab=multiplyAuxsgemmc"sgemmc"ab-- | Matrix product based on BLAS's /cgemm/.multiplyQ::Matrix(ComplexFloat)->Matrix(ComplexFloat)->Matrix(ComplexFloat)multiplyQab=multiplyAuxcgemmc"cgemmc"ab-----------------------------------------------------------------------------foreignimportccallunsafe"svd_l_R"dgesvd::TMMVMforeignimportccallunsafe"svd_l_C"zgesvd::TCMCMVCMforeignimportccallunsafe"svd_l_Rdd"dgesdd::TMMVMforeignimportccallunsafe"svd_l_Cdd"zgesdd::TCMCMVCM-- | Full SVD of a real matrix using LAPACK's /dgesvd/.svdR::MatrixDouble->(MatrixDouble,VectorDouble,MatrixDouble)svdR=svdAuxdgesvd"svdR".fmat-- | Full SVD of a real matrix using LAPACK's /dgesdd/.svdRd::MatrixDouble->(MatrixDouble,VectorDouble,MatrixDouble)svdRd=svdAuxdgesdd"svdRdd".fmat-- | Full SVD of a complex matrix using LAPACK's /zgesvd/.svdC::Matrix(ComplexDouble)->(Matrix(ComplexDouble),VectorDouble,Matrix(ComplexDouble))svdC=svdAuxzgesvd"svdC".fmat-- | Full SVD of a complex matrix using LAPACK's /zgesdd/.svdCd::Matrix(ComplexDouble)->(Matrix(ComplexDouble),VectorDouble,Matrix(ComplexDouble))svdCd=svdAuxzgesdd"svdCdd".fmatsvdAuxfstx=unsafePerformIO$dou<-createMatrixColumnMajorrrs<-createVector(minrc)v<-createMatrixColumnMajorccapp4fmatxmatuvecsmatvstreturn(u,s,transv)wherer=rowsxc=colsx-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'.thinSVDR::MatrixDouble->(MatrixDouble,VectorDouble,MatrixDouble)thinSVDR=thinSVDAuxdgesvd"thinSVDR".fmat-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'.thinSVDC::Matrix(ComplexDouble)->(Matrix(ComplexDouble),VectorDouble,Matrix(ComplexDouble))thinSVDC=thinSVDAuxzgesvd"thinSVDC".fmat-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'.thinSVDRd::MatrixDouble->(MatrixDouble,VectorDouble,MatrixDouble)thinSVDRd=thinSVDAuxdgesdd"thinSVDRdd".fmat-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'.thinSVDCd::Matrix(ComplexDouble)->(Matrix(ComplexDouble),VectorDouble,Matrix(ComplexDouble))thinSVDCd=thinSVDAuxzgesdd"thinSVDCdd".fmatthinSVDAuxfstx=unsafePerformIO$dou<-createMatrixColumnMajorrqs<-createVectorqv<-createMatrixColumnMajorqcapp4fmatxmatuvecsmatvstreturn(u,s,transv)wherer=rowsxc=colsxq=minrc-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'.svR::MatrixDouble->VectorDoublesvR=svAuxdgesvd"svR".fmat-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'.svC::Matrix(ComplexDouble)->VectorDoublesvC=svAuxzgesvd"svC".fmat-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'.svRd::MatrixDouble->VectorDoublesvRd=svAuxdgesdd"svRd".fmat-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'.svCd::Matrix(ComplexDouble)->VectorDoublesvCd=svAuxzgesdd"svCd".fmatsvAuxfstx=unsafePerformIO$dos<-createVectorqapp2gmatxvecsstreturnswherer=rowsxc=colsxq=minrcgracapanbpb=fracapa00nullPtrnbpb00nullPtr-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'.rightSVR::MatrixDouble->(VectorDouble,MatrixDouble)rightSVR=rightSVAuxdgesvd"rightSVR".fmat-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'.rightSVC::Matrix(ComplexDouble)->(VectorDouble,Matrix(ComplexDouble))rightSVC=rightSVAuxzgesvd"rightSVC".fmatrightSVAuxfstx=unsafePerformIO$dos<-createVectorqv<-createMatrixColumnMajorccapp3gmatxvecsmatvstreturn(s,transv)wherer=rowsxc=colsxq=minrcgracapa=fracapa00nullPtr-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'.leftSVR::MatrixDouble->(MatrixDouble,VectorDouble)leftSVR=leftSVAuxdgesvd"leftSVR".fmat-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'.leftSVC::Matrix(ComplexDouble)->(Matrix(ComplexDouble),VectorDouble)leftSVC=leftSVAuxzgesvd"leftSVC".fmatleftSVAuxfstx=unsafePerformIO$dou<-createMatrixColumnMajorrrs<-createVectorqapp3gmatxmatuvecsstreturn(u,s)wherer=rowsxc=colsxq=minrcgracaparucupunbpb=fracaparucupunbpb00nullPtr-----------------------------------------------------------------------------foreignimportccallunsafe"eig_l_R"dgeev::TMMCVMforeignimportccallunsafe"eig_l_C"zgeev::TCMCMCVCMforeignimportccallunsafe"eig_l_S"dsyev::CInt->TMVMforeignimportccallunsafe"eig_l_H"zheev::CInt->TCMVCMeigAuxfstm=unsafePerformIO$dol<-createVectorrv<-createMatrixColumnMajorrrapp3gmatmveclmatvstreturn(l,v)wherer=rowsmgracapa=fracapa00nullPtr-- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/.-- The eigenvectors are the columns of v. The eigenvalues are not sorted.eigC::Matrix(ComplexDouble)->(Vector(ComplexDouble),Matrix(ComplexDouble))eigC=eigAuxzgeev"eigC".fmateigOnlyAuxfstm=unsafePerformIO$dol<-createVectorrapp2gmatmveclstreturnlwherer=rowsmgracapanlpl=fracapa00nullPtrnlpl00nullPtr-- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'.-- The eigenvalues are not sorted.eigOnlyC::Matrix(ComplexDouble)->Vector(ComplexDouble)eigOnlyC=eigOnlyAuxzgeev"eigOnlyC".fmat-- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/.-- The eigenvectors are the columns of v. The eigenvalues are not sorted.eigR::MatrixDouble->(Vector(ComplexDouble),Matrix(ComplexDouble))eigRm=(s',v'')where(s,v)=eigRaux(fmatm)s'=fixeig1sv'=toRows$transvv''=fromColumns$fixeig(toLists')v'eigRaux::MatrixDouble->(Vector(ComplexDouble),MatrixDouble)eigRauxm=unsafePerformIO$dol<-createVectorrv<-createMatrixColumnMajorrrapp3gmatmveclmatv"eigR"return(l,v)wherer=rowsmgracapa=dgeevracapa00nullPtrfixeig1s=toComplex'(subVector0r(asReals),subVectorrr(asReals))wherer=dimsfixeig[]_=[]fixeig[_][v]=[comp'v]fixeig((r1:+i1):(r2:+i2):r)(v1:v2:vs)|r1==r2&&i1==(-i2)=toComplex'(v1,v2):toComplex'(v1,scale(-1)v2):fixeigrvs|otherwise=comp'v1:fixeig((r2:+i2):r)(v2:vs)wherescale=vectorMapValRScalefixeig__=error"fixeig with impossible inputs"-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.-- The eigenvalues are not sorted.eigOnlyR::MatrixDouble->Vector(ComplexDouble)eigOnlyR=fixeig1.eigOnlyAuxdgeev"eigOnlyR".fmat-----------------------------------------------------------------------------eigSHAuxfstm=unsafePerformIO$dol<-createVectorrv<-createMatrixColumnMajorrrapp3fmatmveclmatvstreturn(l,v)wherer=rowsm-- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/.-- The eigenvectors are the columns of v.-- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order).eigS::MatrixDouble->(VectorDouble,MatrixDouble)eigSm=(s',fliprlv)where(s,v)=eigS'(fmatm)s'=fromList.reverse.toList$s-- | 'eigS' in ascending ordereigS'::MatrixDouble->(VectorDouble,MatrixDouble)eigS'=eigSHAux(dsyev1)"eigS'".fmat-- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/.-- The eigenvectors are the columns of v.-- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order).eigH::Matrix(ComplexDouble)->(VectorDouble,Matrix(ComplexDouble))eigHm=(s',fliprlv)where(s,v)=eigH'(fmatm)s'=fromList.reverse.toList$s-- | 'eigH' in ascending ordereigH'::Matrix(ComplexDouble)->(VectorDouble,Matrix(ComplexDouble))eigH'=eigSHAux(zheev1)"eigH'".fmat-- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'.-- The eigenvalues are sorted in descending order.eigOnlyS::MatrixDouble->VectorDoubleeigOnlyS=vrev.fst.eigSHAux(dsyev0)"eigS'".fmat-- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'.-- The eigenvalues are sorted in descending order.eigOnlyH::Matrix(ComplexDouble)->VectorDoubleeigOnlyH=vrev.fst.eigSHAux(zheev1)"eigH'".fmatvrev=flatten.flipud.reshape1-----------------------------------------------------------------------------foreignimportccallunsafe"linearSolveR_l"dgesv::TMMMforeignimportccallunsafe"linearSolveC_l"zgesv::TCMCMCMforeignimportccallunsafe"cholSolveR_l"dpotrs::TMMMforeignimportccallunsafe"cholSolveC_l"zpotrs::TCMCMCMlinearSolveSQAuxfstab|n1==n2&&n1==r=unsafePerformIO$dos<-createMatrixColumnMajorrcapp3fmatamatbmatsstreturns|otherwise=error$st++" of nonsquare matrix"wheren1=rowsan2=colsar=rowsbc=colsb-- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'.linearSolveR::MatrixDouble->MatrixDouble->MatrixDoublelinearSolveRab=linearSolveSQAuxdgesv"linearSolveR"(fmata)(fmatb)-- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'.linearSolveC::Matrix(ComplexDouble)->Matrix(ComplexDouble)->Matrix(ComplexDouble)linearSolveCab=linearSolveSQAuxzgesv"linearSolveC"(fmata)(fmatb)-- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'.cholSolveR::MatrixDouble->MatrixDouble->MatrixDoublecholSolveRab=linearSolveSQAuxdpotrs"cholSolveR"(fmata)(fmatb)-- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'.cholSolveC::Matrix(ComplexDouble)->Matrix(ComplexDouble)->Matrix(ComplexDouble)cholSolveCab=linearSolveSQAuxzpotrs"cholSolveC"(fmata)(fmatb)-----------------------------------------------------------------------------------foreignimportccallunsafe"linearSolveLSR_l"dgels::TMMMforeignimportccallunsafe"linearSolveLSC_l"zgels::TCMCMCMforeignimportccallunsafe"linearSolveSVDR_l"dgelss::Double->TMMMforeignimportccallunsafe"linearSolveSVDC_l"zgelss::Double->TCMCMCMlinearSolveAuxfstab=unsafePerformIO$dor<-createMatrixColumnMajor(maxmn)nrhsapp3fmatamatbmatrstreturnrwherem=rowsan=colsanrhs=colsb-- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'.linearSolveLSR::MatrixDouble->MatrixDouble->MatrixDoublelinearSolveLSRab=subMatrix(0,0)(colsa,colsb)$linearSolveAuxdgels"linearSolverLSR"(fmata)(fmatb)-- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'.linearSolveLSC::Matrix(ComplexDouble)->Matrix(ComplexDouble)->Matrix(ComplexDouble)linearSolveLSCab=subMatrix(0,0)(colsa,colsb)$linearSolveAuxzgels"linearSolveLSC"(fmata)(fmatb)-- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.linearSolveSVDR::MaybeDouble-- ^ rcond->MatrixDouble-- ^ coefficient matrix->MatrixDouble-- ^ right hand sides (as columns)->MatrixDouble-- ^ solution vectors (as columns)linearSolveSVDR(Justrcond)ab=subMatrix(0,0)(colsa,colsb)$linearSolveAux(dgelssrcond)"linearSolveSVDR"(fmata)(fmatb)linearSolveSVDRNothingab=linearSolveSVDR(Just(-1))(fmata)(fmatb)-- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.linearSolveSVDC::MaybeDouble-- ^ rcond->Matrix(ComplexDouble)-- ^ coefficient matrix->Matrix(ComplexDouble)-- ^ right hand sides (as columns)->Matrix(ComplexDouble)-- ^ solution vectors (as columns)linearSolveSVDC(Justrcond)ab=subMatrix(0,0)(colsa,colsb)$linearSolveAux(zgelssrcond)"linearSolveSVDC"(fmata)(fmatb)linearSolveSVDCNothingab=linearSolveSVDC(Just(-1))(fmata)(fmatb)-----------------------------------------------------------------------------------foreignimportccallunsafe"chol_l_H"zpotrf::TCMCMforeignimportccallunsafe"chol_l_S"dpotrf::TMMcholAuxfsta=dor<-createMatrixColumnMajornnapp2fmatamatrstreturnrwheren=rowsa-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.cholH::Matrix(ComplexDouble)->Matrix(ComplexDouble)cholH=unsafePerformIO.cholAuxzpotrf"cholH".fmat-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/.cholS::MatrixDouble->MatrixDoublecholS=unsafePerformIO.cholAuxdpotrf"cholS".fmat-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version).mbCholH::Matrix(ComplexDouble)->Maybe(Matrix(ComplexDouble))mbCholH=unsafePerformIO.mbCatch.cholAuxzpotrf"cholH".fmat-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/ ('Maybe' version).mbCholS::MatrixDouble->Maybe(MatrixDouble)mbCholS=unsafePerformIO.mbCatch.cholAuxdpotrf"cholS".fmat-----------------------------------------------------------------------------------foreignimportccallunsafe"qr_l_R"dgeqr2::TMVMforeignimportccallunsafe"qr_l_C"zgeqr2::TCMCVCM-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/.qrR::MatrixDouble->(MatrixDouble,VectorDouble)qrR=qrAuxdgeqr2"qrR".fmat-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/.qrC::Matrix(ComplexDouble)->(Matrix(ComplexDouble),Vector(ComplexDouble))qrC=qrAuxzgeqr2"qrC".fmatqrAuxfsta=unsafePerformIO$dor<-createMatrixColumnMajormntau<-createVectormnapp3fmatavectaumatrstreturn(r,tau)wherem=rowsan=colsamn=minmn-----------------------------------------------------------------------------------foreignimportccallunsafe"hess_l_R"dgehrd::TMVMforeignimportccallunsafe"hess_l_C"zgehrd::TCMCVCM-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/.hessR::MatrixDouble->(MatrixDouble,VectorDouble)hessR=hessAuxdgehrd"hessR".fmat-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/.hessC::Matrix(ComplexDouble)->(Matrix(ComplexDouble),Vector(ComplexDouble))hessC=hessAuxzgehrd"hessC".fmathessAuxfsta=unsafePerformIO$dor<-createMatrixColumnMajormntau<-createVector(mn-1)app3fmatavectaumatrstreturn(r,tau)wherem=rowsan=colsamn=minmn-----------------------------------------------------------------------------------foreignimportccallunsafe"schur_l_R"dgees::TMMMforeignimportccallunsafe"schur_l_C"zgees::TCMCMCM-- | Schur factorization of a square real matrix, using LAPACK's /dgees/.schurR::MatrixDouble->(MatrixDouble,MatrixDouble)schurR=schurAuxdgees"schurR".fmat-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/.schurC::Matrix(ComplexDouble)->(Matrix(ComplexDouble),Matrix(ComplexDouble))schurC=schurAuxzgees"schurC".fmatschurAuxfsta=unsafePerformIO$dou<-createMatrixColumnMajornns<-createMatrixColumnMajornnapp3fmatamatumatsstreturn(u,s)wheren=rowsa-----------------------------------------------------------------------------------foreignimportccallunsafe"lu_l_R"dgetrf::TMVMforeignimportccallunsafe"lu_l_C"zgetrf::TCMVCM-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/.luR::MatrixDouble->(MatrixDouble,[Int])luR=luAuxdgetrf"luR".fmat-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/.luC::Matrix(ComplexDouble)->(Matrix(ComplexDouble),[Int])luC=luAuxzgetrf"luC".fmatluAuxfsta=unsafePerformIO$dolu<-createMatrixColumnMajornmpiv<-createVector(minnm)app3fmatavecpivmatlustreturn(lu,map(pred.round)(toListpiv))wheren=rowsam=colsa-----------------------------------------------------------------------------------typeTWa=CInt->PD->atypeTQa=CInt->CInt->PC->aforeignimportccallunsafe"luS_l_R"dgetrs::TMVMMforeignimportccallunsafe"luS_l_C"zgetrs::TQ(TW(TQ(TQ(IOCInt))))-- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/.lusR::MatrixDouble->[Int]->MatrixDouble->MatrixDoublelusRapivb=lusAuxdgetrs"lusR"(fmata)piv(fmatb)-- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/.lusC::Matrix(ComplexDouble)->[Int]->Matrix(ComplexDouble)->Matrix(ComplexDouble)lusCapivb=lusAuxzgetrs"lusC"(fmata)piv(fmatb)lusAuxfstapivb|n1==n2&&n2==n=unsafePerformIO$dox<-createMatrixColumnMajornmapp4fmatavecpiv'matbmatxstreturnx|otherwise=error$st++" on LU factorization of nonsquare matrix"wheren1=rowsan2=colsan=rowsbm=colsbpiv'=fromList(map(fromIntegral.succ)piv)::VectorDouble