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

MathLink Example 2: Sum a list of floating point numbers

This example considers a sum of of floating point numbers.
----------------------------------------------------------------
addRealList.c
----------------------------------------------------------------

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



double addRealList( void )
{
 int i,n;

 double *list; 

 MLGetReal64List(stdlink,&list,&n);

 double sum=0;
 for(i=0;i<n;i++) sum = sum + list[i];

 
  MLReleaseReal64List(stdlink,list,n);

 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

------------------------------------------------------
addRealList.tm
------------------------------------------------------

:Begin:
:Function:       addRealList

:Pattern:        addRealList[ L:{___Real}]
:Arguments:      { L  }
:ArgumentTypes:  { Manual  }
:ReturnType:     Real

:End:

:Evaluate: addRealList::usage = "addRealList[x] gives the sum of the elements on x, given that they are real (double floating numbers)"

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

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

addRealList : addRealListtm.c
 ${CXX} addRealListtm.c addRealList.c -o addRealList -lML64i3 -lm -lpthread -lrt -lstdc++

addRealListtm.c : addRealList.tm
 ${MPREP} addRealList.tm -o $@

clean:

 rm addRealListtm.c
 rm addRealList
 

--------------------------------------------------------
addRealList.nb
--------------------------------------------------------

link = Install["./addRealList"]

?addRealList

addRealList[{1., 5., 6.}]

MathLink Example 1: Sum a List of Integers

This example passes a list of integers, performs the sum and returns an integer.
We must note that the length of the list is not type int but type long.

-------------------------------------------------------------------------
addList.c
-------------------------------------------------------------------------

#include "mathlink.h"
extern double addList( double *list, long n);



double addList( double *list, long n)
{

 int i;
 double sum=0;
 for(i=0;i<n;i++) sum = sum + list[i];

 
 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

---------------------------------------------------------------------
 addList.tm
---------------------------------------------------------------------

:Begin:
:Function:       addList

:Pattern:        addList[ L_List]
:Arguments:      { L  }
:ArgumentTypes:  { RealList  }
:ReturnType:     Real

:End:

:Evaluate: addList::usage = "addList[x] gives the sum of the elements of the list x"

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

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

BINARIES = addList

all : $(BINARIES)

addList : addListtm.c addList.c
 ${CXX}  addList.c addListtm.c  -lML64i3 -lm -lpthread -lrt -lstdc++ -o $@

addListtm.c : addList.tm
 ${MPREP} addList.tm -o addListtm.c

-------------------------------------------
addList.nb
-------------------------------------------

link = Install["./addList"]

addList[{1, 5, 6}]

MathLink Example 0

MathLink is a protocol on how to communicate Mathematica with external programs. For example, one would like to run compiled code in C within the Mathematica environment.  I am going to post a sequence of MathLink programs for Linux with an increasing level of complexity, starting from a simple example already found in the Mathematica help (included for completeness) finishing with an example where I send and receive a matrix of complex numbers processed by the GPU using CUDA

Most examples are made of 4 files:

• *.c             : where the computation occurs.

• *.tm          : MathLink template where the input and output are specified and the documentation of the function is placed. 
• Makefile   : where the compilation sequence is specified
• *.nb          : Mathematica notebook that executes the whole program

These files have to be placed in the same folder and compiled with: $make , which may need editing in order to provide the path of the executable file mprep that comes as part of Mathematica. 
  
As seen in Makefile, mprep generates a proper *.c code from the *.tm file, which is compiled along with the *.c file that contains the computational process.

The first example adds two integers

---------------------------------------------------------------------------------
addtwo.c
---------------------------------------------------------------------------------

#include "mathlink.h"
extern int addtwo( int i, int j);



int addtwo( int i, int j)
{

 return i+j;
}



#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

-----------------------------------------------------------------
 addtwo.tm
-----------------------------------------------------------------

int addtwo P(( int, int));

:Begin:
:Function:       addtwo

:Pattern:        AddTwo[i_Integer, j_Integer]
:Arguments:      { i, j }
:ArgumentTypes:  { Integer, Integer }
:ReturnType:     Integer

:End:

:Evaluate: AddTwo::usage = "AddTwo[x, y] gives the sum of two machine integers x and y."

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

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

BINARIES = addtwo

all : $(BINARIES)

addtwo : addtwotm.c addtwo.c
 ${CXX}  addtwotm.c addtwo.c  -lML64i3 -lm -lpthread -lrt -lstdc++ -o $@

addtwotm.c : addtwo.tm
 ${MPREP} addtwo.tm -o addtwotm.c

------------------------------------------------
addtwo.nb
------------------------------------------------

SetDirectory["/home/rcabrera/Documents/source/mathematica/MathLinkExamples/mytest2/alternative"]

link = Install["./addtwo"]

?AddTwo

AddTwo[2, 3]

Subversion

Subversion is a version control software very useful for software developers
Here I found a nice video tutorial
http://showmedo.com/videotutorials/video?name=950000&fromSeriesID=95

Efimov states

I've heard about Efimov states some time ago and now I decided to put them in my shopping list

http://www.sciencedaily.com/releases/2009/12/091211131526.htm

The Challenge to make quantum computers

The challenge to implement a completely scalable  quantum computer is tied with the understanding of the quantum/classical transition where quantum mechanics dominates the explanation of the microscopic world (up to the molecular scale more or less), but classical mechanics explains very well the ordinary macroscopic world. A completely scalable quantum computer would eventually give us a macroscopic taste of the imaginable strange quantum effects that so far have been seen only in the microscopic world.

"Because there are no known fundamental obstacles to such scalability (practical quantum computer with large number of qubits), it has been suggested that failure to achieve it would reveal new physics" -Emanuel Knill


I feel that the most recent papers are becoming more conservative about their predictions on the feasibility of quantum computing despite fact that there is work stating that fault tolerant quantum computation is possible with two ingredients:

• Maximum Error/Gate about 10^-4 to 10^-5
• Effective error correction codes with ancillary qubits.

The last ingredient seems to be more or less accomplished  and the former one does not seem to be fundamentally unattainable despite the fact that we are currently very far. A paper that summarizes these facts with a high degree of scepticism is

Is Fault-Tolerant Quantum Computation Really Possible?

which I also like for its entertaining and straightforward writing style. 

Humble opinion about God and science

For some reason many physicists in the past and present such as Einstein and Hawking used the term God when they wanted to describe something essential and fundamental such as the laws of physics or the Big Bang theory. For most people this was an indication of their religiosity. However, anybody who knows the context of those words can tell that truth is different. Here is an article explaining this issue.

My personal point of view is that science does not require God by definition and not because science can prove that God does not exist. Defining science as the rational understanding of nature based on fundamental laws (and axioms), the acknowledgement of the existence of God would contradict this premise. Using God to explain something in science would be equivalent to admit an irrational element in the theory. In other words, admitting something that is not based on fundamental laws or axioms.

However, I think there is still the inconsistent possibility to put God outside science as the something beyond the laws of physics and axioms. Some would place God as the creator of the laws of physics and axioms. I do not see any consistent and rational way to put God as the explanation of absolutely anything. Having said that, I do not consider myself as completely consistent and rational being and I do not think anybody was, is or will ever be, so, here it lies a window for the endless debate for the existence of God.

Will we ever find a rational and consistent explanation for our consciousness and existence? Something tells me that the most likely answer is NOT.

Inverse Tomography

A few years ago I had the pleasure to attend Pedro Goldman's presentation (Argentinian scientist) about his research in inverse tomography and applications to the calculation of the optimal dosage of radiation for radiotherapy.

Fast Optimization or the Radiation Therapy of Tumors -the Impossible Possible


Sitting down on my chair I would not expect that this was going to be one of the most shocking and interesting presentations I ever attended.

This research field is mathematically interesting and extremely important for medical physics.

Computational Law

I did not know about this this research field until my last trip to Beijing China where I met Set Chandler, who develops his programs in Mathematica

http://www.wolfram.com/products/mathematica/portraits/sethchandler/

The general idea is to develop models as response of the implementation of the laws. The idea is to minimize the subjectivity in today's approach.

A much more ambitious project is to analyse the self consistency of the laws and their level of complexity. The complexity does not only depend on the document's length but also in how intricate is the relation with itself and other laws.

Videos created in Mathematica

I just learned about a video channel for videos created in Mathematica.

http://vimeo.com/channels/mathematica

Medical Imaging

Medical Imaging is today a very active multidisciplinary research field where differential geometry, high performance computing  and diverse mathematics meet together.

In my trip to China I had the opportunity to see some of the latests developments at professor Bart M. ter Haar Romeny's group at the University of Eindhoven in Holland.

Some of the lastest thesis in his group that can be donwloaded are

• E.M. Franken, Enhancement of Crossing Elongated Structures in Images, PhD. Thesis, 2008, Technische Universiteit Eindhoven. Advisors: B.M. ter Haar Romenij, R. Duits

• M.A.van Almsick, Context Models of Lines and Contours, PhD. Thesis, 2007, Technische Universiteit Eindhoven. Advisors: B.M. ter Haar Romeny, L.M.J. Florack

• E. Balmashnova, Scale-Euclidean Invariant Object Retrieval, PhD. Thesis, 2007, Technische Universiteit Eindhoven.Advisors: B.M. ter Haar Romenij,co-advisor: L.M.J. Florack

• F.M.W. Kanters, Towards Object Based Image Editing, PhD. Thesis, 2007, TU/e. Advisors: B.M. ter Haar Romenij, L.M.J. Florack

• B. Platel, Exploring the Deep Structure of Images, PhD. Thesis, 2007, Technische Universiteit Eindhoven. Advisors: B.M. ter Haar Romenij, Co-advisor: L.M.J. Florack

Additionally, professor Romeny has a book entitled


Front-End Vision and Multi-Scale Image Analysis: Multi-scale Computer Vision Theory and Applications, written in Mathematica

Online version

The complete package can be found at
http://extras.springer.com/2003/978-1-4020-1507-6

If you are not afraid of modestly advanced mathematics and want something really applied, try this.

The canonical coset decomposition of unitary matrices through Householder transformations

My paper describing the connection between the Householder decomposition of and the canonical coset parametrization of unitary operators was published today in the Journal of Mathematical Physics

J. Math. Phys. 51, 082101 (2010)

http://arxiv.org/abs/1008.2477

This parametrization is very useful for calculating the parameters that uniquely characterises a given unitary operator U(N). I show that instead of solving a complicated systems of non-linear equations one can perform very stable and efficient Householder transformations.  Right now I am working on some practical applications based on that result.

Independent Component Analysis

I was already familiar with the Principal component analysis, which is a method for finding a linear combination of variables such that the possible correlation between them disappears i.e. elimination of the diagonal elements of the covariance matrix. Statistical independence implies the elimination of correlation but the opposite is not necessarily true because there is possible association between the variables on higher orders.  The Independent Component Analysis seems to be the generalization needed to go beyond the simple correlation.

• Bartlett, M.S. and Movellan, J.R. and Sejnowski, T.J., Face recognition by independent component analysis, IEEE Transactions on neural networks,13, pp 1450--1464, 2002

• Hyvarinen, A. and Oja, E. Independent component analysis: a tutorial, Neural Networks, 13, pp 411--430 2000

10th International Mathematica Symposium

I just returned from the 10th International Mathematica Symposium held in Beijing China, where I gave a talk. It was a great opportunity to see the work of other people around the world and the new features of Mathematica 8.

Some of the topics included medical imaging, stochastic simulations, video editing, real-time control, automatic C code generation, CUDA and more.

Quantum entanglement and classical chaos

I recently attended a talk about the connection between quantum entanglement and classical chaos

http://www.nature.com/nature/journal/v461/n7265/full/nature08396.html

This tells me that the underlying reason of the sensitivity to initial conditions in classical mechanics may be the rise of entanglement in a sequence with the collapse of the wave function.

Venn diagrams in R

This is my first attempt to use R, the well known program for statistics.

I am using the package LIMMA and the instructions found here in order to draw Venn diagrams. My code is

-----------------------------------------
library(limma)

data1<-read.csv("/home/rcabrera/Documents/source/mathematica/Xiongmao/VennDiagram/H1h1.csv", sep=',', header=T)
data2<-read.csv("/home/rcabrera/Documents/source/mathematica/Xiongmao/VennDiagram/H1h2.csv", sep=',', header=T)
data3<-read.csv("/home/rcabrera/Documents/source/mathematica/Xiongmao/VennDiagram/H1h3.csv", sep=',', header=T)

set1<-data1$x
set2<-data2$x
set3<-data3$x

set123<-union(set1,union(set2,set3))

bool1 = logical(length(set123))
bool2 = logical(length(set123))
bool3 = logical(length(set123))

for(i in 1:length(set123)) bool1[i]<-is.element(set123[i],set1)
for(i in 1:length(set123)) bool2[i]<-is.element(set123[i],set2)
for(i in 1:length(set123)) bool3[i]<-is.element(set123[i],set3)

c3<-cbind(bool1, bool2, bool3)

a <- vennCounts(c3)

vennDiagram(a)
-----------------------------------------------

Note.- The data is stored in a column with label "x"

Mersenne Twister in CUDA

The Mersenne Twister algorithm (what a cool name) for generating random numbers is implemented as an example in the SDK folder provided by NVIDIA. However, as it happens with all the examples, it depends on the CUTIL library, which is not part of the standard CUDA API. There is another problem, the compilation of each example is carried out by a general set of makefiles customized for compiling all the examples. These are big problems if we want to generate a stand alone version of the examples. This is desirable because in this way we could adapt the code for our own purposes and eventually link with extra code.

In the rest of this post I will give the makefiles that accomplish the stand alone compilation independent of CUTIL for the Mersenne Twister example.

You will have to follow the following steps :

- Make a copy of the Mersenne Twister folder from the SDK (I renamed as MersenneTwisterCUDA),

- Remove all the references for the CUTIL library

- Add the following lines on MersenneTwister.cu

char raw_path[] = "/YOUR_PATH/MersenneTwisterCUDA/data/MersenneTwister.raw";
char dat_path[] = "/YOUR_PATH/MersenneTwisterCUDA/data/MersenneTwister.dat";

while removing the corresponding lines that originally define raw_path and dat_path.

- Finally, paste the make files I will provide and run $make on the command line
(Define YOUR_PATH according to the structure of your directories)


file: common
----------------------------------------------------------------

CUDA_PATH = /usr/local/cuda


C := gcc
CC := c++
NVCC := $(CUDA_PATH)/bin/nvcc


CUDA_INCLUDE_PATH = /usr/local/cuda/include
CUDA_LIB_PATH := $(CUDA_PATH)/lib64
CUDA_FLAGS = -lcuda -lcudart


===================================================

file: Makefile:
----------------------------------------------------------------------------
include ./common.mk

MersenneTwister.out: genmtrand.c.o MersenneTwister_gold.cpp.o MersenneTwister.cu.o
$(NVCC) -L$(CUDA_LIB_PATH) $(CUDA_FLAGS) genmtrand.c.o MersenneTwister_gold.cpp.o MersenneTwister.cu.o -o MersenneTwister.out

MersenneTwister.cu.o: MersenneTwister.cu dci.h MersenneTwister.h MersenneTwister_kernel.cu
$(NVCC) -c -I$(CUDA_INCLUDE_PATH) MersenneTwister.cu -o MersenneTwister.cu.o

MersenneTwister_gold.cpp.o: MersenneTwister_gold.cpp
$(CC) -c MersenneTwister_gold.cpp -o MersenneTwister_gold.cpp.o

genmtrand.c.o: genmtrand.c
$(C) -c genmtrand.c -o genmtrand.c.o

Classical Origin of Fermion Spin

The spin was postulated by Pauli from experimental evidence, but it was only with the arrival of the Dirac equation that the spin appears naturally. This leaded to many people to consider the spin as fundamentally "quantum". In the following paper we argue that the spin appears naturally from classical relativistic mechanics alone


W. E. Baylis, R Cabrera, J.D. Keselica, Quantum/Classical Interface: Classical Geometric Origin of Fermion Spin, Advances in Applied Clifford Algebras, 2010

Maxwell's demon

Mark G Raizen recently gave a talk about techniques for trapping and cooling atoms. In this way, he actually implemented Maxwell's demon in practice!

http://www.sciencemag.org/cgi/content/abstract/324/5933/1403

Yes, this demon exists and does not violate the second law of thermodynamics because information is a thermodynamic variable with an essential role.

Attention span

The amount of information we can absorb not only depends on the time we actually engage our attention. It also depends on how we distribute this time in subintervals. It seems that we can only pay our highest degree of attention for a few seconds [wikipedia article], so what do we do in the rest?

Related to this topic is the strategy followed by Hollywood, who are interested in maintaining their movies as appealing as possible. In this next article there is an interesting analysis of the duration and distribution of the shoot lengths.

James E. Cutting, Attention and the Evolution of Hollywood Film Psychological Science, 2010

One of the main conclusions is that the distribution of power in the frequency domain must obey 1/f , where f is the frequency. I am sure this is particularly important for teachers and students.

More comments about this article can be found here.

Eigenvalues: C/Lapack

Here is another example on the use of Lapack. This time the objective is to calculate the eigenvectors and eigenvalues of a complex matrix

#include<stdio.h>
#include<math.h>

#include<complex.h>
#include <stdlib.h>

//.......................................................................................................
void zgeTranspose( double complex *Transposed, double complex *M ,int n)
{

int i,j;
for(i=0;i<n;i++)

  for(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}

//......................................................................................................
//  MatrixComplexEigensystem: computes the eigenvectors and eigenValues of input matrix A
//  The eigenvectors are stored in columns
//.....................................................................................................
void MatrixComplexEigensystem( double complex *eigenvectorsVR, double complex *eigenvaluesW, double complex *A, int N)
{

 int i;
  
  double complex *AT = (double complex*) malloc( N*N*sizeof(double complex) );

  zgeTranspose( AT, A , N);
  
  char JOBVL ='N';   // Compute Right eigenvectors

  char JOBVR ='V';   // Do not compute Left eigenvectors

  double complex VL[1];

  int LDVL = 1; 
  int LDVR = N;

 int LWORK = 4*N; 

 double complex *WORK =  (double complex*)malloc( LWORK*sizeof(double complex));

 double complex *RWORK = (double complex*)malloc( 2*N*sizeof(double complex));

 
 int INFO;

 zgeev_( &JOBVL, &JOBVR, &N, AT ,  &N , eigenvaluesW ,

   VL, &LDVL,
   eigenvectorsVR, &LDVR, 
   WORK, 
   &LWORK, RWORK, &INFO );

 zgeTranspose( AT, eigenvectorsVR , N);

 for(i=0;i<N*N;i++) eigenvectorsVR[i]=AT[i];

 
  free(WORK);
  free(RWORK);
  free(AT);
}


int main()
{
  int i,j;
  const int N = 3;

  
  double complex A[] = { 1.+I , 2. ,  3 , 4. , 5.+I , 6. , 7., 8., 9. + I};

  double complex eigenVectors[N*N];
  double complex eigenValues[N];

  
  MatrixComplexEigensystem( eigenVectors, eigenValues, A, N);
  
  printf("\nEigenvectors\n");

  for(i=0;i<N;i++){
    for(j=0;j<N;j++) printf(" (%f,%f) \t", eigenVectors[i*N + j]);

    printf("\n");
  }
  
  printf("\nEigenvalues \n");
  for(i=0;i<N;i++) printf("\n (%f, %f) \t",  eigenValues[i] );

  
printf("\n------------------------------------------------------------\n"); 
return 0;

}

Complex matrix inverse: C++/Lapack

So, here is an example on how to call lapack from c++

#include<iostream>
#include<math.h>

#include<complex>
#include <stdlib.h>

using namespace std;

extern "C" void zgetrf_( int*, int* , complex<double>* , int*, int* , int* );

extern "C" void zgetri_( int*, complex<double>* , int*, int* , complex<double>*, int* , int* );

//........................................................................................
void zgeTranspose( complex<double> *Transposed, complex<double> *M ,int n)
{

int i,j;
for(i=0;i<n;i++)

 for(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}

//.........................................................................................
void MatrixComplexInverse(complex<double> *invA, complex<double> *A, int n)
{

 int LWORK=10*n;

 int *permutations;

 complex<double> *WORK, *tempA;

 tempA = new complex<double>[n*n];

 permutations =  new int[2*n];
 WORK = new complex<double>[n*n];


 int INFO;

 zgeTranspose(tempA,A,n);


 zgetrf_( &n, &n, tempA , &n, permutations , &INFO );

 if (INFO != 0) {
   cout<<"ComplexMatrixInverse: Error at zgetrf  \n"; exit(0);
   }



 zgetri_( &n, tempA , &n, permutations , WORK, &LWORK, &INFO );

 if (INFO != 0) {
   cout<<"ComplexMatrixInverse: Error at zgetri  \n"; exit(0);
   }

 zgeTranspose(invA,tempA,n);

 delete [] WORK;

 delete [] tempA;
 delete [] permutations;

}

/////////////////////////////////////////////////////////////////////////////////////////

int main()
{
 int i,j;
 const int N = 3;

 complex<double> I(0.,1.);

 complex<double> A[] = { 1. + I , 2. ,  3 , 4. , 5.+I , 6. , 7., 8., 9. + I};

 complex<double> invA[N*N];

 MatrixComplexInverse(invA,A,N);


 for(i=0;i<N;i++){
   for(j=0;j<N;j++) cout << invA[i*N + j]<<"\t";

   cout<<"\n";
 }

cout<<"---------------------------\n";

return 0;

}

////////////////////=====================////////////////////

and the make file is

#
CC = c++

#edit LAPACK_PATH if necessary
LAPACK_PATH = /usr/lib64/atlas

a,out: main.cpp
 $(CC)  main.cpp  -L$(LAPACK_PATH) -llapack -lblas  -lgfortran -lm 

My basic rules of programming

Here I am listing some general rules of programming that I learned to follow.
Some of them are enforced by certain languages or programming environments but some are very permissive and require to be more careful. These rules are already well known but here they are again

• Define very well the input from the output arguments of a function. The input arguments should remain unmodified in the routine.
  
• Never write a routine that cannot fit in one page of computer screen in the worst case.
• Develop and validate a program in reasonable small steps. It may take longer to write the code but it will take less time to debug.
• Use a maximum of two nested loops in a certain routine.
• Use meaningful names for the variables even if they become long.
• Do not repeat the same task in different places.
  

Complex matrix inverse: C/Lapack

Here a simple example on how to use Lapack within C to invert matrices using the LU decomposition.

The C++ code is different and I will post that another day

[The html code for this post was generated by http://www.bedaux.net/cpp2html/]

#include<stdio.h>
#include<math.h>

#include<complex.h>
#include <stdlib.h>

//........................................................................................
void zgeTranspose( double complex *Transposed, double complex *M ,int n)
{

int i,j;
for(i=0;i<n;i++)

for(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}

//.........................................................................................
void MatrixComplexInverse(double complex *invA, double complex *A, int n)
{

int LWORK=10*n;

int *permutations;

double complex *WORK, *tempA;

tempA = (double complex*) malloc( n*n*sizeof(double complex) );

permutations = (int*) malloc( 2*n*sizeof(int) );

WORK = (double complex *)malloc(LWORK*sizeof(double complex));


int INFO;

zgeTranspose(tempA,A,n);


zgetrf_( &n, &n, tempA , &n, permutations , &INFO );

if (INFO != 0) {
 printf("ComplexMatrixInverse: Error at zgetrf  \n"); exit(0);
 }



zgetri_( &n, tempA , &n, permutations , WORK, &LWORK, &INFO );

if (INFO != 0) {
 printf("ComplexMatrixInverse: Error at zgetri  \n"); exit(0);
 }

zgeTranspose(invA,tempA,n);

free(WORK);

free(tempA);
free(permutations);

}


/////////////////////////////////////////////////////////////////////////////////////////


int main()
{
int i,j;
const int N = 3;


double complex A[] = { 1.+I , 2. ,  3 , 4. , 5.+I , 6. , 7., 8., 9. + I};

double complex invA[N*N];

MatrixComplexInverse(invA,A,N);


for(i=0;i<N;i++){
 for(j=0;j<N;j++) printf(" (%f,%f) \t", invA[i*N + j]);

 printf("\n");
}

printf("------------------------------------------------------------\n");
return 0;

}


The make file is

#
CC = gcc

#edit LAPACK_PATH if necessary
LAPACK_PATH = /usr/lib64/atlas

a,out: main.c
$(CC)  main.c  -L$(LAPACK_PATH) -llapack -lblas  -lgfortran -lm

Complex matrix inverse III: CUDA

Here is cuinverse_kernel.cu

__global__ void cgeMatrixInverse_kernel(cuFloatComplex *invA , cuFloatComplex *A , int N , cuFloatComplex *Work)
{

int i;
int idx = blockDim.x * blockIdx.x + threadIdx.x;

cuFloatComplex * ThreadWorkSpace;

ThreadWorkSpace = Work + idx*cgeMatrixInverse_WorkSpace()*N*N;

for(i=0;i<N*N;i++) A[ i + idx*N*N ] =  A[i];

A[ idx*N*N ] = make_cuFloatComplex( (float) idx , 1./sqrtf(2));

cgeMatrixInverse(invA + idx*N*N , A + idx*N*N , N , ThreadWorkSpace);


}

_______________________

and my makefile is

#  edit the path of your working directory
BASE_PATH = /home/rcabrera/Documents/source/c/inverseCUDA

#

CUDA_PATH = /usr/local/cuda

CC = gcc
NVCC    := $(CUDA_PATH)/bin/nvcc

CUDA_INCLUDE_PATH = /usr/local/cuda/include
CUDA_LIB_PATH := $(CUDA_PATH)/lib64
CUDA_LIBS = -lcuda -lcudart

a.out: main.cu cuinverse_kernel.cu cumatrixtools.h
$(NVCC)  -I$(CUDA_INCLUDE_PATH) -L$(CUDA_LIB_PATH) $(CUDA_LIBS) -lm  main.cu -o a.out

# $(NVCC) -deviceemu  -I$(CUDA_INCLUDE_PATH) -L$(CUDA_LIB_PATH) $(CUDA_LIBS) -lm  main.cu -o a.out

clean :
rm -f a.out

Complex matrix inverse II: CUDA

So, here is cumatrixtools.h

#include<cuComplex.h>
#ifndef _CUMATRIXTOOLS_H_
#define _CUMATRIXTOOLS_H_


////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ void
cgeSquareMatrixProduct(cuFloatComplex *out, cuFloatComplex *A, cuFloatComplex *B, int N  )
{

int i,j,k;
cuFloatComplex sum;

for(i=0;i<N;i++)

for(j=0;j<N;j++){
  //sum = 0. + 0.*I;
    sum  = make_cuFloatComplex(0.,0.);

  for(k=0;k<N;k++){
    sum =  cuCaddf(sum  , cuCmulf(A[i*N+k],B[k*N+j]) );
  }

  out[i*N+j] = sum;
}


}
//////////////////////////////////////////////////////////////////////

__device__ __host__ static __inline__ void cgeTranspose(cuFloatComplex  *At, cuFloatComplex  *A, int N)
{

int i,j;
for( i = 0; i<N; i++)

for( j = 0; j<N; j++)
  At[i+j*N] = A[i*N+j] ;
}

//////////////////////////////////////////////////////////////////////////////////////
//                LU decomposition of complex matrices
/////////////////////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ void cgeDoolittle_LU_Decomposition(cuFloatComplex *LU, cuFloatComplex *A, int n)
{

 int i, j, k, p;
 cuFloatComplex *p_k, *p_row, *p_col;

for(k=0; k<n*n; k++) LU[k]=A[k];


for (k = 0, p_k = LU; k < n; p_k += n, k++) {

    for (j = k; j < n; j++) {

       for (p = 0, p_col = LU; p < k; p_col += n,  p++)

//           *(p_k + j) -= *(p_k + p) * *(p_col + j);
       *(p_k + j) = cuCsubf( *(p_k + j) , cuCmulf( *(p_k + p) , *(p_col + j) ));
    }

    if ( cuCabsf(*(p_k + k)) != 0.0 )  //return -1;


    for (i = k+1, p_row = p_k + n; i < n; p_row += n, i++) {

       for (p = 0, p_col = LU; p < k; p_col += n, p++)

         // *(p_row + k) -= *(p_row + p) * *(p_col + k);
    *(p_row + k) = cuCsubf( *(p_row + k) , cuCmulf(  *(p_row + p) , *(p_col + k) ));
         //*(p_row + k) /= *(p_k + k);


  *(p_row + k) = cuCdivf( *(p_row + k) , *(p_k + k)  );
    }
 }

}

//////////////////////////////////////////////////////////////////////////////////////////////////
// Back substitution for lower triangular matrices assuming that the diagonal is filled with 1's
// Given  T x = b ,
// where T is a lower triangula matrix NxN with 1's in the diagonal
// b is a known vector
// x is an unknown vector
// the routine otputs x = xout
//////////////////////////////////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ void cgeBackSubstitutionLowerDiagOne(cuFloatComplex *xout, cuFloatComplex *T , cuFloatComplex *b ,int N  )
{

cuFloatComplex bTemp;
int i,j;

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

  bTemp = b[i];
  for(j=0;j<i;j++) bTemp = cuCsubf( bTemp , cuCmulf(T[ i*N + j], xout[j] ) );

  xout[i] = bTemp ;    
}
}

/////////////////////////////////////////////////////////////////////////////////////
// Back substitution for upper triangular matrice
////////////////////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ void cgeBackSubstitutionUpper(cuFloatComplex *xout, cuFloatComplex *T , cuFloatComplex *b ,int N  )
{

cuFloatComplex bTemp;
int i,j;

for(i=N-1;i>=0;i--){

  bTemp = b[i];
  for(j=i+1;j<N;j++) bTemp = cuCsubf( bTemp , cuCmulf(T[ i*N + j], xout[j] ) );

  xout[i] = cuCdivf( bTemp , T[ i*N + i] );    
}
}

///////////////////////
__device__ __host__ static __inline__ void cgeIdentity(cuFloatComplex *Id, int N)
{

int i,j;
for(i=0;i<N;i++)

for(j=0;j<N;j++){
Id[i*N+j] = make_cuFloatComplex(0.f,0.f);
}


for(i=0;i<N;i++){
  Id[i*N+i] = make_cuFloatComplex(1.f,0.f);
}
}

//////////////////////////////////////////////////////////////////////////////////////
//  Inverse of a matrix using the triangularization of matrices
//  Warning:
//          It does not destroys the original matrix A
//          A work space for 3 complex matrices must be supplied in W
//////////////////////////////////////////////////////////////////////////////////////
__device__ __host__ static __inline__ int cgeMatrixInverse_WorkSpace(){ return 3;}

__device__ __host__ static __inline__ void cgeMatrixInverse(cuFloatComplex *inv , cuFloatComplex *A , int N , cuFloatComplex *W)
{

int i;
cuFloatComplex *Id;
cuFloatComplex *invL, *invU;


Id = W;      //Double purpose work space
  invL = W + N*N;

invU = W + 2*N*N;


cgeIdentity( Id, N);

 
 cgeDoolittle_LU_Decomposition(inv,A,N);


for(i=0;i<N;i++)   cgeBackSubstitutionLowerDiagOne( invL + i*N , inv , Id + i*N ,  N  );


for(i=0;i<N;i++) cgeBackSubstitutionUpper( invU + i*N , inv , Id + i*N ,  N  );


cgeSquareMatrixProduct(Id , invL , invU , N  );


cgeTranspose(inv , Id , N);

}


#endif