[comp.lang.fortran] Need complex matrix solver for CM.

john%spectre.unm.edu@ariel.unm.edu (John Prentice) (06/27/91)

I need a routine to solve a dense complex valued linear system on the
CM2.  What I would like is something similiar to what Linpack has, in
terms of functionality, but I will take what I can get.  CMSSL has
code to solve dense real valued linear systems, but not complex valued
ones (please correct me if I am wrong about this).  Any help would be
much appreciated.

John
-- 
John K. Prentice    john@spectre.unm.edu (Internet)
Computational Physics Group
Amparo Corporation, Albuquerque, NM

john@spectre.unm.edu (John Prentice) (06/27/91)

I need a routine to solve a dense complex valued linear system on the
CM2, probably using a LU decomposition and Gaussian elimination.  From
what I know of it, CMSSL only has software like this for real valued
matrices (please correct me if I am wrong), while I need complex valued
ones.  Any help would be most appreciated.

John
-- 
John K. Prentice    john@spectre.unm.edu (Internet)
Computational Physics Group
Amparo Corporation, Albuquerque, NM

francis@hanauma.stanford.edu (Francis Muir) (06/27/91)

John Prentice writes:

  I need a routine to solve a dense complex valued linear system on the
  CM2, probably using a LU decomposition and Gaussian elimination.  From
  what I know of it, CMSSL only has software like this for real valued
  matrices (please correct me if I am wrong), while I need complex valued
  ones.  Any help would be most appreciated.

Re-write the code in REAL; it will probably illuminate.
						
							Fido

francis@hanauma.stanford.edu

john@spectre.unm.edu (John Prentice) (06/28/91)

In article <1991Jun27.161918.10218@leland.Stanford.EDU> francis@hanauma.stanford.edu (Francis Muir) writes:
>John Prentice writes:
>
>  I need a routine to solve a dense complex valued linear system on the
>  CM2, probably using a LU decomposition and Gaussian elimination.  From
>  what I know of it, CMSSL only has software like this for real valued
>  matrices (please correct me if I am wrong), while I need complex valued
>  ones.  Any help would be most appreciated.
>
>Re-write the code in REAL; it will probably illuminate.
>						
>							Fido
>
>francis@hanauma.stanford.edu

This is an amazingly presumptuous, not to mention useless, comment.  If I
have to write it in REAL, I will.  But to do so completely obscures the
elegant structure of the matrices I am solving and requires both artificial
and inelegant constructs.  I would vastly prefer to have software to handle
the complex matrices as they are and keep my code simple.

John
-- 
John K. Prentice    john@spectre.unm.edu (Internet)
Computational Physics Group
Amparo Corporation, Albuquerque, NM

francis@hanauma.stanford.edu (Francis Muir) (06/28/91)

John Prentice writes:

  Francis Muir writes:

    John Prentice writes:

      I need a routine to solve a dense complex valued linear system on the
      CM2, probably using a LU decomposition and Gaussian elimination.  From
      what I know of it, CMSSL only has software like this for real valued
      matrices (please correct me if I am wrong), while I need complex valued
      ones.  Any help would be most appreciated.

    Re-write the code in REAL; it will probably illuminate.
						
  This is an amazingly presumptuous, not to mention useless, comment.  If I
  have to write it in REAL, I will.  But to do so completely obscures the
  elegant structure of the matrices I am solving and requires both artificial
  and inelegant constructs.  I would vastly prefer to have software to handle
  the complex matrices as they are and keep my code simple.

?

If you were to look at:

	CMSSL for CM FORTRAN, Version 2.1, March 1991

you would see in Chapter 2, Linear Algebra Routines, that all routines handle
COMPLEX as well as REAL. For dense systems CMSSL offers both Gauss-Jordan and
QR.

							Fido

bfp@mace.cc.purdue.edu (Bryan Putnam) (06/28/91)

In article <1991Jun27.190724.26161@ariel.unm.edu> john@spectre.unm.edu (John Prentice) writes:
>In article <1991Jun27.161918.10218@leland.Stanford.EDU> francis@hanauma.stanford.edu (Francis Muir) writes:
>>Re-write the code in REAL; it will probably illuminate.
>
>This is an amazingly presumptuous, not to mention useless, comment.  If I

That sounds like a perfectly good solution to me, although it makes the
linear system twice as large. If you need a High School math reminder:

Ax=b ==>

+ A(r)  -A(i) +  x(r)          b(r)
|             |           =    
+ A(i)   A(r) +  x(i)          b(i)

john@spectre.unm.edu (John Prentice) (06/28/91)

In article <7887@mace.cc.purdue.edu> bfp@mace.cc.purdue.edu (Bryan Putnam) writes:
>In article <1991Jun27.190724.26161@ariel.unm.edu> john@spectre.unm.edu (John Prentice) writes:
>>In article <1991Jun27.161918.10218@leland.Stanford.EDU> francis@hanauma.stanford.edu (Francis Muir) writes:
>>>Re-write the code in REAL; it will probably illuminate.
>>
>>This is an amazingly presumptuous, not to mention useless, comment.  If I
>
>That sounds like a perfectly good solution to me, although it makes the
>linear system twice as large. If you need a High School math reminder:
>
>Ax=b ==>
>
>+ A(r)  -A(i) +  x(r)          b(r)
>|             |           =    
>+ A(i)   A(r) +  x(i)          b(i)


Thanks, but I know perfectly well how to convert a complex linear system
to a real linear system.  As I said before, there are reasons I don't wish
to do so having to do with the structure of the matrix I am working with.
If you really need an example of where one might not wish to convert
a complex to a real linear system, consider a diagonally dominant
banded complex matrix.  Look at your real matrix above.  Both the 
A(r) and A(i) block matrices will now be diagonally dominant and banded,
but the big matrix will not be.  By converting this to a real matrix,
you lose both the diagonal domiance and you introduce off diagonal bands.  
Both make solving the system alot harder with anything but Gaussian 
elimination.  My systems are order 10,000 to 100,000.  You can't solve
them with straight Gaussian elimination and expect to get the answers back
anytime soon.  That means you need iterative schemes and these often
depend on the structure of the matrix such as the bandedness and diagonal
dominance.  One could construct schemes that would work with the
real version of the matrix, but why do it if software exists to do the problem
as is?  It just makes the problem obscure.  In my case, I am solving my 
system using an iterative scheme which at some point exploits the fact 
that my matrices have diagonal bands which are large compared to the 
rest of the matrix.  I solve the smaller banded problem using Gaussian 
elimination as part of my iterative scheme.  Reworking the algorithms and 
proving their convergence properties is not worth the effort if complex 
matrix solvers exist, as they in fact appear to do.

What bugged me about the original comment was not the suggestion to convert
to a real matrix, it was the comment "it will probably illuminate".  In
either case, my question has been answered, so perhaps we can let this
particular discussion end.

John
-- 
John K. Prentice    john@spectre.unm.edu (Internet)
Computational Physics Group
Amparo Corporation, Albuquerque, NM

winstead@faraday.ece.cmu.edu (Charles Holden Winstead) (06/28/91)

>In article <1991Jun27.190724.26161@ariel.unm.edu> john@spectre.unm.edu (John Prentice) writes:
>linear system twice as large. If you need a High School math reminder:
                                             ^^^^^^^^^^^^^^^^^^^^^^^^^^

Yo dude, grow up!  I agree that it's a simple task of changing a complex 
linear system into a real one, but with a language that deals with complex
numbers as nicely as fortran does, let's agree that there are some benefits
to not switching from complex to real and back again repeatedly.  Why not
just look through the code that works for real matrices and change the 
corresponding variable declaration from real to complex.

>
>Ax=b ==>
>
>+ A(r)  -A(i) +  x(r)          b(r)
>|             |           =    
>+ A(i)   A(r) +  x(i)          b(i)

Gee thanks.

-Chuck

john@spectre.unm.edu (John Prentice) (06/29/91)

In article <1991Jun28.061454.29461@fs7.ece.cmu.edu> winstead@faraday.ece.cmu.edu (Charles Holden Winstead) writes:
>>In article <1991Jun27.190724.26161@ariel.unm.edu> john@spectre.unm.edu (John Prentice) writes:
>>linear system twice as large. If you need a High School math reminder:
>                                             ^^^^^^^^^^^^^^^^^^^^^^^^^^
>
>Yo dude, grow up!  I agree that it's a simple task of changing a complex 
>linear system into a real one, but with a language that deals with complex
>numbers as nicely as fortran does, let's agree that there are some benefits
>to not switching from complex to real and back again repeatedly.  Why not
>just look through the code that works for real matrices and change the 
>corresponding variable declaration from real to complex.
>

The first time you come across a piece of code in a matrix solver that
checks to see what sign a real number is, how are you going to change this
to handle complex numbers?  It takes more than just changing the variable
declarations.  In either case, it is a sterile argument.  I judge it to
be preferential in my situation to use complex numbers and I don't 
particularly feel compelled to argue about it.  I am certainly not 
proselytizing that this is the only way to do it.

Peace!

John
-- 
John K. Prentice    john@spectre.unm.edu (Internet)
Computational Physics Group
Amparo Corporation, Albuquerque, NM

bfp@mace.cc.purdue.edu (Bryan Putnam) (06/29/91)

In article <1991Jun28.061454.29461@fs7.ece.cmu.edu> winstead@faraday.ece.cmu.edu (Charles Holden Winstead) writes:
>
>Yo dude, grow up!  I agree that it's a simple task of changing a complex 
>linear system into a real one, but with a language that deals with complex
>numbers as nicely as fortran does, let's agree that there are some benefits
>to not switching from complex to real and back again repeatedly.  Why not
>

There can also be some benefits to separating the data into real and
imaginary parts, because on some machines the vector hardware
instructions do not operate on complex data. If you the programmer do
not separate the data, then the compiler has too.  For example, the
highly optimized FFT routines written for the CYBER 205, and which
were used by the National Weather Service, start off by SEPARATING
the real and imaginary parts simply because the vector hardware is
more efficient at performing real arithmetic. Too many programmers
get into trouble by ignoring the architecture of specific machines.

>Why not
>just look through the code that works for real matrices and change the 
>corresponding variable declaration from real to complex.

I suppose the obvious answer would be that the source code is seldom
available, unless of course this is a trick question.

jbs@WATSON.IBM.COM (06/29/91)

         Rewriting a dense complex system of linear equations as a real
system will not only use twice the storage it will take twice as long to
solve (using Gaussian elimination).  It is comparable to solving a sym-
metric system with a general system solver.
                          James B. Shearer