Lapack to MPL
Well I've spent at least 2 months on f2c trying to modify it so that the output is c++ using MpIeee instead of double. We got pretty close and I learned a lot but we seem to be stuck with parameter passing between functions (inside the f2c compiler the type information of parameters is not always available and it is very cryptic/oldstyle c code so modification of f2c is tedious).
Today while doing my usual sourceforge project browsing I stumbled on f2matlab. Which could turn out to be a key component in getting LAPACK to mpl. I'm writing this down not to forget it by tomorrow since it is all speculary but it's a very interesting idea, here's a little schematic:
LAPACK (f77 code)
-> to_f90 tool
-> LAPACK (f90 code)
-> f2matlab (expects f90 files)
-> m-files already usable in MPL but we do some more steps
-> MPL compiler
-> generates C++ files using our own multi-precision matrices based on boost
-> create library
-> use what we've learned in the precompiler GSL api-generator to
generate code for use inside compiler and interpreter of MPL
After these steps we should now have working builtin multi-precision
kernel routines generated from lapack !
What are the obstacles, first of al to_f90 by Alan Miller. It requires a proprietary compiler for fortran90. I'm trying out the intel fortran compiler tomorrow, just to get the linux executable from to_f90.f90. Next is hoping that f2matlab delivers. It looks very promising but we have to wait and see if it's capable of converting all lapack routines.
Next step is getting f2matlab (which is oddly enough a compiler written in matlab itself) to run in MPL, getting the whole system almost self hosting (only the to_f90 remains a external component). In the end if all goes well we'll be able to package all these steps into an automated system that just takes in a fortran file (f77 or f90) and spits out a C++ compiled binary module for use in MPL. That would be so sweat... Ok time to go home and have dinner, can't wait for tomorrow :D (5/11/2006, and I'm starving because it's 20h30pm).
News update today (12/05/2006), I was able to convert these files already to .m files:
wschrep@pascal:~/fortran2mpl/lapack-mpieee/matlab-lapack/lapack> ls
dasum.m dgesv.m dlagtf.m dlarzb.m dorgrq.m dsbev.m dsyrfs.m
daxpy.m dgesvx.m dlagtm.m dlarz.m dorgtr.m dsbevx.m dsyrk.m
dbdsdc.m dgetc2.m dlagts.m dlarzt.m dorm2l.m dsbgst.m dsyr.m
dbdsqr.m dgetf2.m dlagv2.m dlas2.m dorm2r.m dsbgvd.m dsysv.m
dcopy.m dgetrf.m dlahqr.m dlascl.m dormbr.m dsbgv.m dsysvx.m
ddisna.m dgetri.m dlahrd.m dlasd0.m dormhr.m dsbgvx.m dsytd2.m
ddot.m dgetrs.m dlaic1.m dlasd1.m dorml2.m dsbmv.m dsytf2.m
dgbbrd.m dggbak.m dlaln2.m dlasd2.m dormlq.m dsbtrd.m dsytrd.m
dgbcon.m dggbal.m dlals0.m dlasd3.m dormql.m dscal.m dsytrf.m
dgbequ.m dgges.m dlalsa.m dlasd4.m dormqr.m dsecnd.m dsytri.m
dgbmv.m dggesx.m dlalsd.m dlasd5.m dormr2.m dspcon.m dsytrs.m
dgbrfs.m dggev.m dlamch.m dlasd6.m dormr3.m dspevd.m dtbcon.m
dgbsv.m dggevx.m dlamrg.m dlasd7.m dormrq.m dspev.m dtbmv.m
dgbsvx.m dggglm.m dlangb.m dlasd8.m dormrz.m dspevx.m dtbrfs.m
dgbtf2.m dgghrd.m dlange.m dlasd9.m dormtr.m dspgst.m dtbsv.m
dgbtrf.m dgglse.m dlangt.m dlasda.m dpbcon.m dspgvd.m dtbtrs.m
dgbtrs.m dggqrf.m dlanhs.m dlasdq.m dpbequ.m dspgv.m dtgevc.m
dgebak.m dggrqf.m dlansb.m dlasdt.m dpbrfs.m dspgvx.m dtgex2.m
dgebal.m dggsvd.m dlansp.m dlaset.m dpbstf.m dspmv.m dtgexc.m
dgebd2.m dggsvp.m dlanst.m dlasq1.m dpbsv.m dspr2.m dtgsen.m
dgebrd.m dgtcon.m dlansy.m dlasq2.m dpbsvx.m dsprfs.m dtgsja.m
dgecon.m dgtrfs.m dlantb.m dlasq3.m dpbtf2.m dspr.m dtgsna.m
dgeequ.m dgtsv.m dlantp.m dlasq4.m dpbtrf.m dspsv.m dtgsy2.m
dgees.m dgtsvx.m dlantr.m dlasq5.m dpbtrs.m dspsvx.m dtgsyl.m
dgeesx.m dgttrf.m dlanv2.m dlasq6.m dpocon.m dsptrd.m dtpcon.m
dgeev.m dgttrs.m dlapll.m dlasr.m dpoequ.m dsptrf.m dtpmv.m
dgeevx.m dgtts2.m dlapmt.m dlasrt.m dporfs.m dsptri.m dtprfs.m
dgegs.m dhgeqz.m dlapy2.m dlassq.m dposv.m dsptrs.m dtpsv.m
dgegv.m dhsein.m dlapy3.m dlasv2.m dposvx.m dstebz.m dtptri.m
dgehd2.m dhseqr.m dlaqgb.m dlaswp.m dpotf2.m dstedc.m dtptrs.m
dgehrd.m dlabad.m dlaqge.m dlasy2.m dpotrf.m dstegr.m dtrcon.m
dgelq2.m dlabrd.m dlaqp2.m dlasyf.m dpotri.m dstein.m dtrevc.m
dgelqf.m dlacon.m dlaqps.m dlatbs.m dpotrs.m dsteqr.m dtrexc.m
dgelsd.m dlacpy.m dlaqsb.m dlatdf.m dppcon.m dsterf.m dtrmm.m
dgels.m dladiv.m dlaqsp.m dlatps.m dppequ.m dstevd.m dtrmv.m
dgelss.m dlae2.m dlaqsy.m dlatrd.m dpprfs.m dstev.m dtrrfs.m
dgelsx.m dlaebz.m dlaqtr.m dlatrs.m dppsv.m dstevr.m dtrsen.m
dgelsy.m dlaed0.m dlar1v.m dlatrz.m dppsvx.m dstevx.m dtrsm.m
dgemm.m dlaed1.m dlar2v.m dlatzm.m dpptrf.m dswap.m dtrsna.m
dgemv.m dlaed2.m dlarfb.m dlauu2.m dpptri.m dsycon.m dtrsv.m
dgeql2.m dlaed3.m dlarfg.m dlauum.m dpptrs.m dsyevd.m dtrsyl.m
dgeqlf.m dlaed4.m dlarf.m dnrm2.m dptcon.m dsyev.m dtrti2.m
dgeqp3.m dlaed5.m dlarft.m dopgtr.m dpteqr.m dsyevr.m dtrtri.m
dgeqpf.m dlaed6.m dlarfx.m dopmtr.m dptrfs.m dsyevx.m dtrtrs.m
dgeqr2.m dlaed7.m dlargv.m dorg2l.m dptsv.m dsygs2.m dtzrqf.m
dgeqrf.m dlaed8.m dlarnv.m dorg2r.m dptsvx.m dsygst.m dtzrzf.m
dgerfs.m dlaed9.m dlarrb.m dorgbr.m dpttrf.m dsygvd.m idamax.m
dger.m dlaeda.m dlarre.m dorghr.m dpttrs.m dsygv.m ieeeck.m
dgerq2.m dlaein.m dlarrf.m dorgl2.m dptts2.m dsygvx.m lsame.m
dgerqf.m dlaev2.m dlarrv.m dorglq.m drotg.m dsymm.m lsamen.m
dgesc2.m dlaexc.m dlartg.m dorgql.m drot.m dsymv.m xerbla.m
dgesdd.m dlag2.m dlartv.m dorgqr.m drscl.m dsyr2k.m
dgesvd.m dlags2.m dlaruv.m dorgr2.m dsbevd.m dsyr2.m
I mailed the author of f2matlab (Benjamin Barrowes) the remaining 4 lapack files that don't convert yet:
wschrep@pascal:~/fortran2mpl/lapack-mpieee> ls bugs
dlacon.f dlaln2.f dlaruv.f ilaenv.f
dlacon.f90 dlaln2.f90 dlaruv.f90 ilaenv.f90
Hope he replies, otherwise will have to convert these files manually which is not such a big deal thus so far progress is looking good... ok time to go home it's 20h pm again, have dinner and skate a bit because it's weekend :)
Just some usefull links for me:
Matlab 7 notes and reference to Benjamin Barrowes.
...
15/5/2006:
After studying the converted f90 code by f2matlab I've notised some problems which need further investigation.
For now there are some obstacles to overcome using this method, the largest seems to be fortran 'go to' statements. Since matlab does not have these, f2matlab does not convert them correctly. The labels are lost in conversion, not the goto statements themselves, but putting the label in some formatted command line like '%L10 CONTINUE' there is a workaround possible. By extending mpl to support goto's during the conversion to c++ we can convert the routines and later disable the goto statement once the conversion is done. Other smaller problems include varargin addition to the parameter list even when this does not seem necessary and external functions are disabled by listing them as local variables (uncommenting these luckily solves it, but this is ofcourse a bit of manual work or another translation pass has to be done after f2matlab).
Mailed author of f2matlab with some questions regarding the implementation.
Some numbers:
wschrep@pascal:~/fortran2mpl/lapack-mpieee/matlab-lapack/lapack> grep "go to" * | cut -d ":" -f -1 | uniq | wc -l
139
wschrep@pascal:~/fortran2mpl/lapack-mpieee/matlab-lapack/lapack> grep "go to" * | cut -d ":" -f -1 | wc -l
587
This means 139 fortran files contain 587 goto statements. A bit much to check manually but if I have to it is possible to do this. Some easy goto's are the ones used to implement a while statement:
Some fortran code :
!+ WHILE( C.EQ.ONE )LOOP
10 CONTINUE
IF( c == one ) THEN
a = 2*a
c = dlamc3( a, one )
c = dlamc3( c, -a )
GO TO 10
END IF
!+ END WHILE
matlab code after f2matlab pass:
%+ WHILE(C.EQ.ONE)LOOP
% continue;
if(c == one) ;
a = 2.*a;
c = dlamc3(a, one);
c = dlamc3(c, -a);
go to 10;
end;
%+ END WHILE
Code after manual change:
while(c == one) ;
a = 2.*a;
c = dlamc3(a, one);
c = dlamc3(c, -a);
end;
But thinking it is always this easy, is a bit optimistic. One good example of tedious goto's is the ilaenv.f fortran file (which was also one of the files that did not convert at all using f2matlab). It contains a switch/case with goto's to various parts in the file, changing this manually for 100 files is too much work...
So what's next, well it depends on Benjamin Barrowes reaction to the mail. But mostly I'm thinking about either a rewrite/modification of his translator by myself or having a last attempt at getting f2c modification finished (again, feel free to look at the f2c sources, you'll quickly agree it's not a trivial job, to modify this legacy/ancient c-code...)
Yet another possible solution is the plusFORT tool, which claims to eliminate many goto's in fortran code: plusFort intro
2006/05/17:
Got positive mail back from Benjamin regarding f2matlab, he fixed most bugs and even did my feature request (to leave the labels in comments so that I can use the goto's).
But in parallel I'm also considering a modification of the f2j translator. This is written in a modern programming style and might be easier to modify for MPL purposes.
Either way, one of these tools will get the job done :).
2006/06/03 Work on a new/improved matrix rep is going well. But we have some serious setbacks with the f2matlab tool. I mailed the author with my problem, and I hope there will be a fix for it:
Hi,
I have stumbled on some more serious conversion problems with f2matlab.
The gotos work without problems now but conversion of subroutine CALL's
don't work as expected.
For example in dgetrf90.f:
CALL dgemm( 'No transpose', 'No transpose', m-j-jb+1, &
n-j-jb+1, jb, -one, a( j+jb, j ), lda, &
a( j, j+jb ), lda, one, a( j+jb, j+jb ), lda )
We call the dgemm subroutine. But f2matlab copies all input args to output
arguments and this gives problems here with expressions being passed as input args:
[dumvar1,dumvar2,dumvar3,dumvar4, jb,dumvar6, a(j+jb, j), lda,a(j, j+jb),
lda, one, a(j+jb, j+jb), lda]
=dgemm('No transpose', 'No transpose', fix(m)-fix(j)-fix(jb)+1,
fix(n)-fix(j)-fix(jb)+1, fix(jb), -one, a(fix(j)+fix(jb), fix(j)),fix(lda),
a(fix(j), fix(j)+fix(jb)), fix(lda), one, a(fix(j)+fix(jb), fix(j)+fix(jb)), fix(lda));
expressions such as a(j+jb,j) cannot be on left side of assignment obviously.
Do you know of any workarounds for this problem?
In the attachment are the sample files (f90 and converted dgetrf.m file).
At first glance a fix seems to be removing the expressions given as input
arguments to functions out of the list of the output
arguments (since these can't be changed/returned by
the called function/sub anyway).
Kind regards,
Walter