Calling Matlab functions in Fortran code - need help
Maku Chan
My work is on numerical techniques which require to solve sparse matrices etc. All my work is in fortran and I have almost no skill in Matlab.But I have just realized that Matlab can do some computations quite easily with its built-in functions and I want to take advantage of this. I have no time to learn Matlab deeply or work on matlab-fortran communication (very much off-focus issues and have deadline soon)
Some details:
I need to solve Ax=b system where A is very large and sparse matrix (20000x20000 or more). At every time step (iteration) I need to solve this system, get 'x' and reuse it for further computation in that step. [A] is fixed and formed before the iteration begins - but vector 'b' changes at every iteration. So in the iteration -- i) vector 'b' is formed, ii) solve Ax=b to find 'x' iii) use many other (explicit algebraic) equations to calculate some other values / fields (magnetic fields) iv) save field values for that step and move to the next step. Then these steps are repeated once and again. I am doing steps i), iii), iv) in fortran and want to do step ii) in matlab. Actually I have the whole code written in fortran but want to do some analysis by just doing the step ii) in Matlab. Can you please advise me if that is possible and if, how can I do that?
I have tried Matlab to solve Ax=b just for a single iteration. That means, I could read matrix [A], an instance of vector 'b' in matlab (after saving them in files from fortran). Then I solved them in Matlab command window using 'bicgstab' command (also used 'luinc' preconditioner etc, did some matrix blocking etc). But I need to repeatedly do this for the whole system and also calculate 'b' vector in fortran using the obtained solution. Can you please advise me on this. Thank you so much.
James Tursa
>
> I need to solve Ax=b system where A is very large and sparse matrix (20000x20000 or more). At every time step (iteration) I need to solve this system, get 'x' and reuse it for further computation in that step. [A] is fixed and formed before the iteration begins - but vector 'b' changes at every iteration.
You could use the MATLAB engine for this. Look in the doc for the External Interfaces API section. You would be using engOpen, engPutVariable, engEvalString, and engGetVariable. But you have to build the mxArray inside your Fortran code first, and mxArray sparse matrix building is not trivial. How are you building the sparse matrix in your Fortran code now?
James Tursa
Maku Chan
Thank you for the reply.
>
> You could use the MATLAB engine for this. Look in the doc for the External Interfaces API section. You would be using engOpen, engPutVariable, engEvalString, and engGetVariable. But you have to build the mxArray inside your Fortran code first, and mxArray sparse matrix building is not trivial. How are you building the sparse matrix in your Fortran code now?
>
Feeling intimidated !
In fortran the sparse matrix comes from a huge number of long equations. I have implemented all those long equations in fortran which produces a large sparse matrix (A) and does all the steps I mentioned before. (thats why it is not possible to re-implement it in matlab)
As my Matlab skill is almost nil I have no idea what you meant by
Matlab engine, and these commands - engOpen, engPutVariable, engEvalString, and engGetVariable. May be I need to look into the documentation.
It sounds not so straightforward. But as I have time constraint I need to keep focussed and just want to take advantage of matlab, if any. I cannot afford to spend long time to learn matlab for this. Do you know any easiest way. Thank you very much.
James Tursa
The engOpen, engPutVariable, engEvalString, and engGetVariable are not very difficult to use. I could easily put together a short demo, or you could look at demos in the MATLAB doc, or you could download some examples from the FEX. The biggest hurdle is the sparse matrix itself. You mention that your Fortran code already produces a sparse matrix. How do you know it is sparse? Is the data stored in some special sparse way? If so, what are the details about this storage? You cannot hope to interface with MATLAB and use MATLAB functions until you know these answers. So there is no point in starting this project until you know how the Fortran data is stored.
James Tursa
Maku Chan
>
From the equations I know it is sparse and I produced the sparsity pattern in fortran (using some available routine). When I read the sparse matrix file in Matlab I could see the same sparsity pattern using 'spy' command. There is no doubt about it as I am working with sparse matrix in fortran everyday.
My sparse matrix data in fortran is available in COO (coordinate), CSR (compressed sparse row) and CCS (compressed column sparse) formats. In fact I convert from one format to another as required by the solver or preconditioner routines (in fortran). Thank you very much.
James Tursa
>
> From the equations I know it is sparse and I produced the sparsity pattern in fortran (using some available routine). When I read the sparse matrix file in Matlab I could see the same sparsity pattern using 'spy' command. There is no doubt about it as I am working with sparse matrix in fortran everyday.
OK. Can you post a very brief Fortran example of a sparse matrix that is in the format that MATLAB accepts? It sounds like you just call some pre-canned routines to do the format conversions for you ... I don't really care about these details. I am just interested in example source code that shows the variables used to store the sparse matrix that you know works with MATLAB when you write it out to a file and then read it into MATLAB. Once that is known, writing the engine code to call MATLAB functions will not be too hard.
James Tursa
John Densmore
I would have to think that using imsl would be faster than matlab.
Maku Chan
I was told matlab function 'load' can read sparse matrix from a file. But I actually read a large nxn sparse matrix in matlab in a different way. I first stored the sparse matrix in COO format (in a file) which stores like the following:
RowIndex ColInDex Value
2 3 0.00001
23 5 0.003033
20 10 0.23333
12 9 0.022
That means all the elements except the above are zeros for this 25x25 matrix. Then I added the header 25 25 4 (m , n , number of nonzeros) in this file. When the file is prepared like above I can read it in Matlab using the following code.
http://math.nist.gov/MatrixMarket/src/rdcoord.m
>> a=rdcoord('a-coo-file');
You can read any sparse matrix in matlab this way as long as the file is properly prepared (which I can easily do from my fortran code). For the vector 'b' I first stored it using fortran code and then read it in Matlab using following commands:
>> format long eng
>> fid=fopen('b-vector-file');
>> b=fscanf(fid,'%e');
>> fclose all
Now we have both 'a' and 'b' in matlab. But, as I wrote before, 'b' changes at every iteration and get prepared in fortran using the results from Ax=b. So the above reading and solving Ax=b need to be done in matlab iteratively again and again.
>(you wrote)
MATLAB when you write it out to a file and then read it into MATLAB. Once that is known, writing the engine code to call MATLAB functions will not be too hard.
>
Now do you think the major difficulty, as you said, is gone? Is it now possible to prepare the code so that one doesnot need to care about the Matlab side -- rather only concerns with the fortran code?
Thank you so much
Bruno Luong
> James,
> I was told matlab function 'load' can read sparse matrix from a file. But I actually read a large nxn sparse matrix in matlab in a different way. I first stored the sparse matrix in COO format (in a file) which stores like the following:
> RowIndex ColInDex Value
> 2 3 0.00001
> 23 5 0.003033
> 20 10 0.23333
> 12 9 0.022
> That means all the elements except the above are zeros for this 25x25 matrix. Then I added the header 25 25 4 (m , n , number of nonzeros) in this file. When the file is prepared like above I can read it in Matlab using the following code.
>
> http://math.nist.gov/MatrixMarket/src/rdcoord.m
>
> >> a=rdcoord('a-coo-file');
>
> You can read any sparse matrix in matlab this way as long as the file is properly prepared (which I can easily do from my fortran code). For the vector 'b' I first stored it using fortran code and then read it in Matlab using following commands:
>
> >> format long eng
> >> fid=fopen('b-vector-file');
> >> b=fscanf(fid,'%e');
> >> fclose all
>
> Now we have both 'a' and 'b' in matlab. But, as I wrote before, 'b' changes at every iteration and get prepared in fortran using the results from Ax=b. So the above reading and solving Ax=b need to be done in matlab iteratively again and again.
Indeed, if you accept interfacing through ascii file (slow and few digits are gone), then things get easier. One you read the matrix, use SPARSE command to build Matlab sparse matrix. Depending on the property of your matrix, you might want factorize with cholesky, lu or qr. This files from File Exchange will do the job for you.
http://www.mathworks.com/matlabcentral/fileexchange/24119
Then solve the linear equation for every iteration.
This should take a couple of lines of code.
Bruno