Inverse Matrix by the Exchange Method

After some time off, here I am again.

Here is an interesting implementation of the inverse of a matrix that I found useful when I was developing a CUDA program. This is the most efficient method to implement the inverse of a matrix in terms of memory usage, which is handy when we want to put everything in the fast but limited shared memory.

The Mathematica code is here

InverseExchange.cdf

The cuda code is next
InverseExchange.cu

J. Math. Phys Most Downloaded Articles in August 2010

My paper about the connection between the Householder decomposition of unitary matrices and the canonical coset representation of unitary groups was the second most downloaded paper in August 2010

http://jmp.aip.org/features/most_downloaded?month=8&year=2010

Ito Calculus with Tensorial

Here is a short example on how to use the Tensorial package for Ito calculus (stochastic calculus)

ito.pdf

Numerical Linear Algebra and Quantum Control

Here is a poster about some of my recently published research

Poster pdf

MathLink example 4: Trace of a complex matrix

In this example, I introduce the ability to treat complex numbers.

-------------------------------------------------------
myTrace.c
-------------------------------------------------------

#include<stdio.h>
#include<string.h>

#include "mathlink.h"
extern void myTrace( void );

void myTrace( void )
{

 int i,j,n;
  
 double *matrix;

 int *dims;
 char **heads;
 int rank; 

 if( !MLGetReal64Array(stdlink,&matrix,&dims,&heads,&rank) ){ 
 // unable to read data

  printf(" MLGetReal64Array error reading data \n");
 };

 n = dims[0];

 double resum=0.;
 double imsum=0.;

 double sum[2];

 int outdims[1];

 char **outheads;
 int outrank = 1; 

 
 for(i=0;i<n;i++){

 resum = resum + matrix[2*i*n + 2*i];

 imsum = imsum + matrix[2*i*n + 2*i + 1];
 }

 
 MLReleaseReal64Array(stdlink,matrix,dims,heads,rank);
 
 MLPutFunction(stdlink,"Complex",2);

 MLPutFloat(stdlink,resum);
 MLPutFloat(stdlink,imsum);


}


#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif

-----------------------------------------------------
myTrace.tm
-----------------------------------------------------

:Begin:
:Function:       myTrace

:Pattern:        myTrace[ L_List]
:Arguments:      { L  }
:ArgumentTypes:  { Manual  }
:ReturnType:     Manual

:End:

:Evaluate: myTrace::usage = "myTrace[M] gives the trace of complex matrix M"

------------------------------------------------------
Makefile
------------------------------------------------------

MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

myTrace : myTracetm.c myTrace.c
 ${CXX} myTracetm.c myTrace.c -o myTrace -lML64i3 -lm -lpthread -lrt -lstdc++

myTracetm.c : myTrace.tm
 ${MPREP} myTrace.tm -o $@

---------------------------------------------------
myTrace.nb
---------------------------------------------------

link = Install["./myTrace"]

(m = {{1., 2., 1.}, {2., 2., 2.}, {3., 3., 3.}} + I) // MatrixForm

myTrace[m]

MathLink example 3: Trace of a real matrix

In this example we calculate the trace of a matrix of floating point numbers (double type in C)

---------------------------------------------------
myTrace.c
---------------------------------------------------

#include "mathlink.h"
extern double myTrace( void );

double myTrace( void )
{
 int i,j,n;
  
 double *matrix;
 int *dims;
 char **heads;

 int rank; 

 if( !MLGetReal64Array(stdlink,&matrix,&dims,&heads,&rank) ){ 
 // unable to read data

  return 0.;
 };

 n = dims[0];

 double sum=0;
 for(i=0;i<n;i++) sum = sum + matrix[i*n+i];
 
 MLReleaseReal64Array(stdlink,matrix,dims,heads,rank);

 return sum;
}

#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif

----------------------------------------------
myTrace.tm
----------------------------------------------

:Begin:
:Function:       myTrace

:Pattern:        myTrace[ L_List]
:Arguments:      { L  }
:ArgumentTypes:  { Manual  }
:ReturnType:     Real

:End:

:Evaluate: myTrace::usage = "myTrace[M] gives the trace of matrix M"

---------------------------------------------
Makefile
---------------------------------------------

MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

myTrace : myTracetm.c myTrace.c
 ${CXX} myTracetm.c myTrace.c -o myTrace -lML64i3 -lm -lpthread -lrt -lstdc++

myTracetm.c : myTrace.tm
 ${MPREP} myTrace.tm -o $@

Catenary

The catenary is the curve that is naturally formed when a chain of constant linear density hangs from two fixed points. This is a typical example on variational calculus showing the use of Lagrange multipliers. Next, I develop the problem entirely in Mathematica

Catenary
Catenary(pdf)

You may need MathematicaPlayer if you do not have Mathematica.

This problem can be generalized to find the surface with minimum energy, which is shown next