Re: 2nd Try!
- To: mathgroup at smc.vnet.net
- Subject: [mg9793] Re: 2nd Try!
- From: Daniel Lichtblau <danl>
- Date: Fri, 28 Nov 1997 05:35:08 -0500
- Organization: Wolfram Research, Inc.
- Sender: owner-wri-mathgroup at wolfram.com
hashemi at sun.iust.ac.ir wrote:
>
> Dear Mr. "Seanross":
>
> Thank you very much for your respond. However there are some more
> questions in which I greatly appreciate your futher advice.
>
> 1) Unfortunatly my matrices are too large (at least 500*500 all complex
> elements!) Furthermore their elements are all products of Bessel
> Functions
> of complex argument; Hence I am not sure if rewriting all my (2000
> lines)
> Fortran program in mathematica will even run on even a 586 pentium?!
> In addition in my Fortran program I invert these large complex
> matrices
> using full pivoting technique with quadraple precision! So I am not
> so
> sure that mathematica can do such large matrix inversion!?
> Please correct me on above, if I am wrong.
>
> 2) Do you know about the program INTERCALL which is suppose to interface
> mathematica with Fortran programs? If so can you please send me
> information
> on how to get access to this software?
>
> Again I sincerely thank you time and efforts and will be very happy to
> hear from you soon. Please remmember "like everthingelse here, our
> email system sometimes does not work properly" Hence I appreciate your
> patience.
>
> Thank you
> Dr. S. M. Hasheminejad
> Mechanical Engineering Dept.
> Iran University of Science & Tech
> Narmak, Tehran 16844 Iran
Here are examples of numeric matrix inversion in Mathematica using
100x100 matrices of random complex values in the rectangle from -10-10I
to 10+10I. First I show timings for the bignum case (35 digits
precision), then the machine arithmetic case. Note that the result is
quite good (it is really only when the matrix is ill-conditioned that
we expect otherwise).
(* bignum case *)
In[4]:= mat = Table[Random[Complex, {-10-10I,10+10I}, 35], {100},
{100}];
In[5]:= Timing[inv=Inverse[mat];]
Out[5]= {244.66 Second, Null}
In[6]:= id=mat.inv];
In[7]:= Max[Abs[id-IdentityMatrix[100]]]
-37
Out[7]= 0. 10
(* machine number case *)
In[9]:= mat2 = Table[Random[Complex, {-10-10I,10+10I}], {100}, {100}];
In[10]:= Timing[inv2=Inverse[mat2];] Out[10]= {0.55 Second, Null}
In[11]:= Timing[id2=mat2.inv2;]
Out[11]= {0.94 Second, Null}
In[17]:= Max[Abs[id2-IdentityMatrix[100]]]
-14
Out[17]= 2.91837 10
Some remarks:
The inversion is performed internally using Gaussian elimination with
partial pivoting. Mathematica uses the standard O(n^3) algorithm (where
n = dimension of matrix), hence you can estimate timings for 500x500
case. For the machine number case I illustrate below.
In[20]:= mat3 = Table[Random[Complex, {-10-10I,10+10I}], {400}, {400}];
In[21]:= Timing[inv3=Inverse[mat3];] Out[21]= {38.44 Second, Null}
We expected it to take around 4^3==64 times longer, and indeed it did.
One caveat is that as size increases you can expect to lose more time
due to effects such as disk I/O.
For most realistic purposes one does not want to invert a matrix, but
rather to solve linear systems of equations. The relevant functions
would be LinearSolve and LUDecomposition/LUBackSubstitution. These are
likely to be numerically more stable as well as faster than using
Inverse and Dot to achieve the same purpose.
Daniel Lichtblau
Wolfram Research