Matrix inversion on MacOS using Accelerate framework

wilbert · Post by **wilbert** » Tue Jun 26, 2018 6:34 am

ImportC "-framework Accelerate"
  dgetrf_(*m, *n, *a, *lda, *ipiv, *info)
  dgetri_(*m, *a, *lda, *ipiv, *work, *lwork, *info)
  sgetrf_(*m, *n, *a, *lda, *ipiv, *info)
  sgetri_(*m, *a, *lda, *ipiv, *work, *lwork, *info)
EndImport


Procedure.l MatrixInvertF(N, *Matrix.Float)
  Protected Error.l, *Pivot, *Workspace 

  *Pivot = AllocateMemory(N<<2, #PB_Memory_NoClear)
  *Workspace = AllocateMemory(N<<2, #PB_Memory_NoClear)
  
  sgetrf_(@N, @N, *Matrix, @N, *Pivot, @Error)
  If Error = 0
    sgetri_(@N, *Matrix, @N, *Pivot, *Workspace, @N, @Error)
  EndIf
  
  FreeMemory(*Pivot)
  FreeMemory(*Workspace)
  
  ProcedureReturn Error
EndProcedure


Procedure.l MatrixInvertD(N, *Matrix.Double)
  Protected Error.l, *Pivot, *Workspace 
  
  *Pivot = AllocateMemory(N<<2, #PB_Memory_NoClear)
  *Workspace = AllocateMemory(N<<3, #PB_Memory_NoClear)
  
  dgetrf_(@N, @N, *Matrix, @N, *Pivot, @Error)
  If Error = 0
    dgetri_(@N, *Matrix, @N, *Pivot, *Workspace, @N, @Error)
  EndIf
  
  FreeMemory(*Pivot)
  FreeMemory(*Workspace)
  
  ProcedureReturn Error
EndProcedure




Dim A.f(3,3)
A(0,0) = 2
A(0,1) = 5
A(0,2) = 3
A(0,3) = 2

A(1,0) = 4
A(1,1) = 10
A(1,2) = 1
A(1,3) = 7

A(2,0) = 1
A(2,1) = 5
A(2,2) = 2
A(2,3) = 2

A(3,0) = 2
A(3,1) = 1
A(3,2) = 2
A(3,3) = 1

MatrixInvertF(4, @A())

For t=0 To 3
  For v=0 To 3
    Debug StrF(A(t,v),4)
  Next
Next