From: Ian Scott <ian.scott@st...>  20030513 10:11:51

> Hi All: > > I am not sure if I have fallen victim to your antispam effort, I am > reposting my query. Since my last post, I have obtained and compiled > clapack. At this point I must decide between using the old VisSDK > library and integrating clapack and an mpeg library, or > integrating the > required routines into vxl's vnl library. I suspect, you have fallen victim to "nobody around on Monday who really knows the answer." We wrapped the fortran public maths libraries in vnl_algo for a reason  most of us don't want to worry about it! BTW  You may as well keep this discussion to vxlmaintainers. I'll have a go at your questions, but I'm not one of the original authors, so I may not be of much help. We are always interested in having useful new routines added to VXL, so I'll try and get some more help. > >  Original Message  > > Subject: Integration of portion of LAPACK/CLAPACK routines to vnl > > Date: Fri, 09 May 2003 15:39:08 0400 > > From: Jeff Strickrott <jstric01@...> > > To: > vxlusers@...,vxlmaintainers@... > > > > Hello All: > > > > I am new to the vxl environment and I am trying to port > code that was > > written to use the NAG numerical library and/or a subset of > LAPACK and > > need to have the following LAPACK/CLAPACK or NAG routines: > > > >  DSYEV  compute all eigenvalues and, optionally, eigenvectors of a > > real symmetric matrix A. > >  DGETRF  compute an LU factorization of a general MbyN matrix A > > using partial pivoting with row interchanges. > >  DGETRI  compute the inverse of a matrix using the LU > factorization > > computed by DGETRF. > >  DGEEV  compute for an NbyN real nonsymmetric matrix A, the > > eigenvalues and, optionally, the left and/or right eigenvectors. I > > believe that your vnl_real_eigensystem will do the same thing. > > > >  F04AEF: calculates the accurate solution of a set of real linear > > equations with multiple righthand sides using an LU > factorization with > > partial pivoting, and iterative refinement. That is given a > set of real > > linear equations AX = B , the routine first computes an LU > factorization > > of A with partial pivoting, PA =LU , where P is a > permutation matrix, L > > is lower triangular and U is unit upper triangular. An > approximation to > > X is found by forward and backward substitution. The > residual matrix R > > =B  AX is then calculated using additional precision, and > a correction > > D to X is found by solving LUD = PR . X is replaced by X + > D and this > > iterative refinement of the solution is repeated until full machine > > accuracy has been obtained. > > > > My expertise is not with computer vision or numerical > routines (yet :) > > ) and I have the following questions: > > 1. Would it be easier to just obtain CLAPACK and use the appropriate > > routines from there? Has anyone compiled both CLAPACK and > vnl together? I'm pretty sure that the CLAPACK library is a trivial syntacticsugar interface to LAPACK. It replaces some of the byreference parameters with byvalue ones. When I did some experiments (Comparing the speed of VNL's matrixvector multiplications and the stuff in Intel's optimised BLAS library MKL (which BTW, gcc+vnl won)http://sourceforge.net/mailarchive/message.php?msg_id=3628964) I found it easier to use the ordinary f2c interfaces  which are identical to the original fortran ones. > > 2. Add the FORTRAN routines per the procedures in the appendix? It would be great to add the fortran routines, and write an vnl_algo interface. You could model the C++ interface design for the LU decomposition on the others, e.g. vnl_svd and vnl_cholesky. > > 3. What are (if any) the thread safety issues with vxl, and in > > particular the vnl libraries? I'm pretty sure that LINPACK is not thread safe. In general threadsafety is an "aim" of vxl, but you may find some code (like vnl_algo/netlib) that isn't. All of the code in vnl (the classes like vnl_vector, etc. but not the stuff in vnl_algo) is thread safe if you disable vnl's own memory manager. I can't see any problems adding any other potentially nonthread safe routines to netlib. Although if it isn't too hard making them thread safe (in the sense that two threads can call the same routine at the same time) that would be great. > > 5. Any other suggestions? To give you some explanation, most of the complex routines in vnl are imports of the older LINPACK. I think the authors found it easier to use the LINPACK routines  since they had less dependencies. Replacing all the LINPACK routines with full LAPACK ones is something that needs to be done  but which noone has yet had the time to do. I am not aware of any problems with LINPACK and LAPACK routines coexisting. They share the the same BLAS routines, but as I mentioned LAPACK needs more helper routines than the authors converted for the existing LINPACK code. Have a look at how vnl/algo/vnl_real_eigensystem.cxx to see how vnl accesses the contents of the netlib library. One other point  it might be worth asking if you really nead the LU decomposition. Would SVD be suitable instead? Finally  If you do want to add the routines to VNL, get yourself a user account on sourceforge, and let me know. I'll set up write access to the repository for you. Ian. 