{-# LANGUAGE ForeignFunctionInterface #-}{-# LANGUAGE FlexibleContexts #-}{-# LANGUAGE FlexibleInstances #-}{-# LANGUAGE BangPatterns #-}------------------------------------------------------------------------------- |-- Module : Data.Packed.Internal.Matrix-- Copyright : (c) Alberto Ruiz 2007-- License : GPL-style---- Maintainer : Alberto Ruiz <aruiz@um.es>-- Stability : provisional-- Portability : portable (uses FFI)---- Internal matrix representation--------------------------------------------------------------------------------- #hidemoduleData.Packed.Internal.Matrix(Matrix(..),rows,cols,cdat,fdat,MatrixOrder(..),orderOf,createMatrix,mat,cmat,fmat,toLists,flatten,reshape,Element(..),trans,fromRows,toRows,fromColumns,toColumns,matrixFromVector,subMatrix,liftMatrix,liftMatrix2,(@@>),atM',saveMatrix,singleton,size,shSize,conformVs,conformMs,conformVTo,conformMTo)whereimportData.Packed.Internal.CommonimportData.Packed.Internal.SignaturesimportData.Packed.Internal.VectorimportForeign.Marshal.Alloc(alloca,free)importForeign.Marshal.Array(newArray)importForeign.Ptr(Ptr,castPtr)importForeign.Storable(Storable,peekElemOff,pokeElemOff,poke,sizeOf)importData.Complex(Complex)importForeign.C.TypesimportForeign.C.String(newCString)importSystem.IO.Unsafe(unsafePerformIO)-----------------------------------------------------------------{- Design considerations for the Matrix Type
-----------------------------------------
- we must easily handle both row major and column major order,
for bindings to LAPACK and GSL/C
- we'd like to simplify redundant matrix transposes:
- Some of them arise from the order requirements of some functions
- some functions (matrix product) admit transposed arguments
- maybe we don't really need this kind of simplification:
- more complex code
- some computational overhead
- only appreciable gain in code with a lot of redundant transpositions
and cheap matrix computations
- we could carry both the matrix and its (lazily computed) transpose.
This may save some transpositions, but it is necessary to keep track of the
data which is actually computed to be used by functions like the matrix product
which admit both orders.
- but if we need the transposed data and it is not in the structure, we must make
sure that we touch the same foreignptr that is used in the computation.
- a reasonable solution is using two constructors for a matrix. Transposition just
"flips" the constructor. Actual data transposition is not done if followed by a
matrix product or another transpose.
-}dataMatrixOrder=RowMajor|ColumnMajorderiving(Show,Eq)transOrderRowMajor=ColumnMajortransOrderColumnMajor=RowMajor{- | Matrix representation suitable for GSL and LAPACK computations.
The elements are stored in a continuous memory array.
-}dataMatrixt=Matrix{irows::{-# UNPACK #-}!Int,icols::{-# UNPACK #-}!Int,xdat::{-# UNPACK #-}!(Vectort),order::!MatrixOrder}-- RowMajor: preferred by C, fdat may require a transposition-- ColumnMajor: preferred by LAPACK, cdat may require a transpositioncdat=xdatfdat=xdatrows::Matrixt->Introws=irowscols::Matrixt->Intcols=icolsorderOf::Matrixt->MatrixOrderorderOf=order-- | Matrix transpose.trans::Matrixt->MatrixttransMatrix{irows=r,icols=c,xdat=d,order=o}=Matrix{irows=c,icols=r,xdat=d,order=transOrdero}cmat::(Elementt)=>Matrixt->Matrixtcmatm@Matrix{order=RowMajor}=mcmatMatrix{irows=r,icols=c,xdat=d,order=ColumnMajor}=Matrix{irows=r,icols=c,xdat=transdatardc,order=RowMajor}fmat::(Elementt)=>Matrixt->Matrixtfmatm@Matrix{order=ColumnMajor}=mfmatMatrix{irows=r,icols=c,xdat=d,order=RowMajor}=Matrix{irows=r,icols=c,xdat=transdatacdr,order=ColumnMajor}-- C-Haskell matrix adapter-- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) rmat::(Storablet)=>Matrixt->(((CInt->CInt->Ptrt->t1)->t1)->IOb)->IObmataf=unsafeWith(xdata)$\p->doletmg=dog(fi(rowsa))(fi(colsa))pfm{- | Creates a vector by concatenation of rows
@\> flatten ('ident' 3)
9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@
-}flatten::Elementt=>Matrixt->Vectortflatten=xdat.cmattypeMtts=Int->Int->Ptrt->s-- not yet admitted by my haddock version-- infixr 6 ::>-- type t ::> s = Mt t s-- | the inverse of 'Data.Packed.Matrix.fromLists'toLists::(Elementt)=>Matrixt->[[t]]toListsm=splitEvery(colsm).toList.flatten$m-- | Create a matrix from a list of vectors.-- All vectors must have the same dimension,-- or dimension 1, which is are automatically expanded.fromRows::Elementt=>[Vectort]->MatrixtfromRowsvs=casecompatdim(mapdimvs)ofNothing->error"fromRows applied to [] or to vectors with different sizes"Justc->reshapec.join.map(adaptc)$vswhereadaptcv|dimv==c=v|otherwise=constantD(v@>0)c-- | extracts the rows of a matrix as a list of vectorstoRows::Elementt=>Matrixt->[Vectort]toRowsm=toRows'0wherev=flattenmr=rowsmc=colsmtoRows'k|k==r*c=[]|otherwise=subVectorkcv:toRows'(k+c)-- | Creates a matrix from a list of vectors, as columnsfromColumns::Elementt=>[Vectort]->MatrixtfromColumnsm=trans.fromRows$m-- | Creates a list of vectors from the columns of a matrixtoColumns::Elementt=>Matrixt->[Vectort]toColumnsm=toRows.trans$m-- | Reads a matrix position.(@@>)::Storablet=>Matrixt->(Int,Int)->tinfixl9@@>m@Matrix{irows=r,icols=c}@@>(i,j)|safe=ifi<0||i>=r||j<0||j>=cthenerror"matrix indexing out of range"elseatM'mij|otherwise=atM'mij{-# INLINE (@@>) #-}-- Unsafe matrix access without range checkingatM'Matrix{icols=c,xdat=v,order=RowMajor}ij=v`at'`(i*c+j)atM'Matrix{irows=r,xdat=v,order=ColumnMajor}ij=v`at'`(j*r+i){-# INLINE atM' #-}------------------------------------------------------------------matrixFromVectorocv=Matrix{irows=r,icols=c,xdat=v,order=o}where(d,m)=dimv`quotRem`cr|m==0=d|otherwise=error"matrixFromVector"-- allocates memory for a new matrixcreateMatrix::(Storablea)=>MatrixOrder->Int->Int->IO(Matrixa)createMatrixordrc=dop<-createVector(r*c)return(matrixFromVectorordcp){- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = trans . reshape r@
where r is the desired number of rows.)
@\> reshape 4 ('fromList' [1..12])
(3><4)
[ 1.0, 2.0, 3.0, 4.0
, 5.0, 6.0, 7.0, 8.0
, 9.0, 10.0, 11.0, 12.0 ]@
-}reshape::Storablet=>Int->Vectort->Matrixtreshapecv=matrixFromVectorRowMajorcvsingletonx=reshape1(fromList[x])-- | application of a vector function on the flattened matrix elementsliftMatrix::(Storablea,Storableb)=>(Vectora->Vectorb)->Matrixa->MatrixbliftMatrixfMatrix{icols=c,xdat=d,order=o}=matrixFromVectoroc(fd)-- | application of a vector function on the flattened matrices elementsliftMatrix2::(Elementt,Elementa,Elementb)=>(Vectora->Vectorb->Vectort)->Matrixa->Matrixb->MatrixtliftMatrix2fm1m2|not(compatm1m2)=error"nonconformant matrices in liftMatrix2"|otherwise=caseorderOfm1ofRowMajor->matrixFromVectorRowMajor(colsm1)(f(xdatm1)(flattenm2))ColumnMajor->matrixFromVectorColumnMajor(colsm1)(f(xdatm1)((xdat.fmat)m2))compat::Matrixa->Matrixb->Boolcompatm1m2=rowsm1==rowsm2&&colsm1==colsm2------------------------------------------------------------------{- | Supported matrix elements.
This class provides optimized internal
operations for selected element types.
It provides unoptimised defaults for any 'Storable' type,
so you can create instances simply as:
@instance Element Foo@.
-}class(Storablea)=>ElementawheresubMatrixD::(Int,Int)-- ^ (r0,c0) starting position ->(Int,Int)-- ^ (rt,ct) dimensions of submatrix->Matrixa->MatrixasubMatrixD=subMatrix'transdata::Int->Vectora->Int->Vectoratransdata=transdataP-- transdata'constantD::a->Int->VectoraconstantD=constantP-- constant'instanceElementFloatwheretransdata=transdataAuxctransFconstantD=constantAuxcconstantFinstanceElementDoublewheretransdata=transdataAuxctransRconstantD=constantAuxcconstantRinstanceElement(ComplexFloat)wheretransdata=transdataAuxctransQconstantD=constantAuxcconstantQinstanceElement(ComplexDouble)wheretransdata=transdataAuxctransCconstantD=constantAuxcconstantC-------------------------------------------------------------------transdata'::Storablea=>Int->Vectora->Int->Vectoratransdata'c1vc2=ifnoneedthenvelseunsafePerformIO$dow<-createVector(r2*c2)unsafeWithv$\p->unsafeWithw$\q->doletgo(-1)_=return()go!i(-1)=go(i-1)(c1-1)go!i!j=dox<-peekElemOffp(i*c1+j)pokeElemOffq(j*c2+i)xgoi(j-1)go(r1-1)(c1-1)returnwwherer1=dimv`div`c1r2=dimv`div`c2noneed=r1==1||c1==1-- {-# SPECIALIZE transdata' :: Int -> Vector Double -> Int -> Vector Double #-}-- {-# SPECIALIZE transdata' :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double) #-}-- I don't know how to specialize...-- The above pragmas only seem to work on top level defs-- Fortunately everything seems to work using the above class-- C versions, still a little faster:transdataAuxfunc1dc2=ifnoneedthendelseunsafePerformIO$dov<-createVector(dimd)unsafeWithd$\pd->unsafeWithv$\pv->fun(fir1)(fic1)pd(fir2)(fic2)pv//check"transdataAux"returnvwherer1=dimd`div`c1r2=dimd`div`c2noneed=r1==1||c1==1transdataP::Storablea=>Int->Vectora->Int->VectoratransdataPc1dc2=ifnoneedthendelseunsafePerformIO$dov<-createVector(dimd)unsafeWithd$\pd->unsafeWithv$\pv->ctransP(fir1)(fic1)(castPtrpd)(fisz)(fir2)(fic2)(castPtrpv)(fisz)//check"transdataP"returnvwherer1=dimd`div`c1r2=dimd`div`c2sz=sizeOf(d@>0)noneed=r1==1||c1==1foreignimportccallunsafe"transF"ctransF::TFMFMforeignimportccallunsafe"transR"ctransR::TMMforeignimportccallunsafe"transQ"ctransQ::TQMQMforeignimportccallunsafe"transC"ctransC::TCMCMforeignimportccallunsafe"transP"ctransP::CInt->CInt->Ptr()->CInt->CInt->CInt->Ptr()->CInt->IOCInt----------------------------------------------------------------------constant'vn=unsafePerformIO$dow<-createVectornunsafeWithw$\p->doletgo(-1)=return()go!k=pokeElemOffpkv>>go(k-1)go(n-1)returnw-- C versionsconstantAuxfunxn=unsafePerformIO$dov<-createVectornpx<-newArray[x]app1(funpx)vecv"constantAux"freepxreturnvconstantF::Float->Int->VectorFloatconstantF=constantAuxcconstantFforeignimportccallunsafe"constantF"cconstantF::PtrFloat->TFconstantR::Double->Int->VectorDoubleconstantR=constantAuxcconstantRforeignimportccallunsafe"constantR"cconstantR::PtrDouble->TVconstantQ::ComplexFloat->Int->Vector(ComplexFloat)constantQ=constantAuxcconstantQforeignimportccallunsafe"constantQ"cconstantQ::Ptr(ComplexFloat)->TQVconstantC::ComplexDouble->Int->Vector(ComplexDouble)constantC=constantAuxcconstantCforeignimportccallunsafe"constantC"cconstantC::Ptr(ComplexDouble)->TCVconstantP::Storablea=>a->Int->VectoraconstantPan=unsafePerformIO$doletsz=sizeOfav<-createVectornunsafeWithv$\p->doalloca$\k->dopokekacconstantP(castPtrk)(fin)(castPtrp)(fisz)//check"constantP"returnvforeignimportccallunsafe"constantP"cconstantP::Ptr()->CInt->Ptr()->CInt->IOCInt------------------------------------------------------------------------ | Extracts a submatrix from a matrix.subMatrix::Elementa=>(Int,Int)-- ^ (r0,c0) starting position ->(Int,Int)-- ^ (rt,ct) dimensions of submatrix->Matrixa-- ^ input matrix->Matrixa-- ^ resultsubMatrix(r0,c0)(rt,ct)m|0<=r0&&0<rt&&r0+rt<=(rowsm)&&0<=c0&&0<ct&&c0+ct<=(colsm)=subMatrixD(r0,c0)(rt,ct)m|otherwise=error$"wrong subMatrix "++show((r0,c0),(rt,ct))++" of "++show(rowsm)++"x"++show(colsm)subMatrix''(r0,c0)(rt,ct)cv=unsafePerformIO$dow<-createVector(rt*ct)unsafeWithv$\p->unsafeWithw$\q->doletgo(-1)_=return()go!i(-1)=go(i-1)(ct-1)go!i!j=dox<-peekElemOffp((i+r0)*c+j+c0)pokeElemOffq(i*ct+j)xgoi(j-1)go(rt-1)(ct-1)returnwsubMatrix'(r0,c0)(rt,ct)(Matrix{icols=c,xdat=v,order=RowMajor})=Matrixrtct(subMatrix''(r0,c0)(rt,ct)cv)RowMajorsubMatrix'(r0,c0)(rt,ct)m=trans$subMatrix'(c0,r0)(ct,rt)(transm)---------------------------------------------------------------------------- | Saves a matrix as 2D ASCII table.saveMatrix::FilePath->String-- ^ format (%f, %g, %e)->MatrixDouble->IO()saveMatrixfilenamefmtm=docharname<-newCStringfilenamecharfmt<-newCStringfmtleto=iforderOfm==RowMajorthen1else0app1(matrix_fprintfcharnamecharfmto)matm"matrix_fprintf"freecharnamefreecharfmtforeignimportccallunsafe"matrix_fprintf"matrix_fprintf::PtrCChar->PtrCChar->CInt->TM----------------------------------------------------------------------conformMsms=map(conformMTo(r,c))mswherer=maximum(maprowsms)c=maximum(mapcolsms)conformVsvs=map(conformVTon)vswheren=maximum(mapdimvs)conformMTo(r,c)m|sizem==(r,c)=m|sizem==(1,1)=reshapec(constantD(m@@>(0,0))(r*c))|sizem==(r,1)=repColscm|sizem==(1,c)=repRowsrm|otherwise=error$"matrix "++shSizem++" cannot be expanded to ("++showr++"><"++showc++")"conformVTonv|dimv==n=v|dimv==1=constantD(v@>0)n|otherwise=error$"vector of dim="++show(dimv)++" cannot be expanded to dim="++shownrepRowsnx=fromRows(replicaten(flattenx))repColsnx=fromColumns(replicaten(flattenx))sizem=(rowsm,colsm)shSizem="("++show(rowsm)++"><"++show(colsm)++")"