Fortran - Who uses clapack?

Author

Who uses clapack?

Victor Eijkhout

2004-12-10, 3:58 pm

The authors of lapack/scalapack are starting to work on a new release of
these packages. One of the things we want to address is the sillyness
that C users, instead of linking to the binaries of the Fortran version,
use an automatically translated C version.

Therefore we'd like to know precisely what the reasons for this are.

Naively one would say, compile the Fortran version, append an underscore
(when appropriate) to routine names, and just link those libraries to
your object files.

What are the gotchas, and what are the real stumbling blocks here?
Please spell it out in as much detail as you can muster. This situation
really needs to improve.

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

beliavsky@aol.com

2004-12-10, 3:58 pm

see.sig@for.addy (Victor Eijkhout) wrote:
>The authors of lapack/scalapack are starting to work on a new release of
>these packages. One of the things we want to address is the sillyness
>that C users, instead of linking to the binaries of the Fortran version,
>use an automatically translated C version.
>
>Therefore we'd like to know precisely what the reasons for this are.
>
>Naively one would say, compile the Fortran version, append an underscore
>(when appropriate) to routine names, and just link those libraries to
>your object files.
>
>What are the gotchas, and what are the real stumbling blocks here?
>Please spell it out in as much detail as you can muster. This situation
>really needs to improve.

The reasons may be more political than technical. At my company, the IT department
is Fortran-phobic, so if someone in the research department hands them a
code with mixed C++ and Fortran 77, they will want the F77 part to be f2c'ed.
They believe that FORTRAN is unmaintainable because they don't know FORTRAN.

I can use Fortran to get work done, but only if I stay out of their sight.
I hate IT people.

I know of at least one other company in my industry (financial services)
with a policy that all production code must be in C++. Researchers can experiment
in Matlab etc., but eventually code must be migrated.

I find it odd that you think people ought to use Fortran LAPACK rather than
clapack, when it appears from http://www.fortranstatement.com/cgi-bin/signatures.pl
that you signed a "petition to retire Fortran", which contains lies such
as "FORTRAN2003 Adds No Functionality".

----== Posted via Newsfeeds.Com - Unlimited-Uncensored-Secure Usenet News==----
http://www.newsfeeds.com The #1 Newsgroup Service in the World! >100,000 Newsgroups
---= 19 East/West-Coast Specialized Servers - Total Privacy via Encryption =---

Tino

2004-12-10, 3:58 pm

Victor Eijkhout

2004-12-10, 3:58 pm

"beliavsky@aol.com" <beliavsky@127.0.0.1:7501> wrote:

> I find it odd that you think people ought to use Fortran LAPACK rather
> than clapack, when it appears from
> http://www.fortranstatement.com/cgi-bin/signatures.pl that you signed a
> "petition to retire Fortran", which contains lies such as "FORTRAN2003
> Adds No Functionality".

Hee. I forgot about that. And I don't think it's relevant for my
question.

Lapack/Scalapack is 500,000 lines of tremendously useful legacy code,
and it's not going to go away anytime soon. At the moment my only
ambition is to make that code as cleanly usable from other languages as
possible. Using automatically translated versions is at the least a
support nightmare.

Thanks for your response,

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Victor Eijkhout

2004-12-10, 3:58 pm

Tino <tino52@yahoo.com> wrote:

> On a slightly related note, is there any effort (or even thought) being
> put into the development of a "standardized" C++ object-oriented LAPACK
> library? Something to wrap around the BLAS and LAPACK, perhaps?

The thought has occurred to many people. I mean, passing LDA into a C++
routine is the programming equivalent of turning a crank handle on the
space shuttle. Now, what to do about it....

I've started looking at Sidl/Babel to get a more native looking
interface in other languages than F77. Not sure to what extent that will
satisfy you.

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Dan Pop

2004-12-10, 3:58 pm

In <1gokq2z.1qtze21fbcwaoN%see.sig@for.addy> see.sig@for.addy (Victor Eijkhout) writes:

>The authors of lapack/scalapack are starting to work on a new release of
>these packages. One of the things we want to address is the sillyness
>that C users, instead of linking to the binaries of the Fortran version,
>use an automatically translated C version.
>
>Therefore we'd like to know precisely what the reasons for this are.
>
>Naively one would say, compile the Fortran version, append an underscore
>(when appropriate) to routine names, and just link those libraries to
>your object files.
>
>What are the gotchas, and what are the real stumbling blocks here?
>Please spell it out in as much detail as you can muster. This situation
>really needs to improve.

The issue that springs to mind is that Fortran and C store matrices
differently: Fortran does it column-wise (the first index varies fastest),
C does it row-wise (the last index varies fastest).

The fact that Fortran array indexing is 1-based (by default) and C array
indexing is 0-based shouldn't cause any problems when using the compiled
Fortran code from C.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de
Currently looking for a job in the European Union

Gordon Sande

2004-12-10, 3:58 pm

Victor Eijkhout wrote:
> "beliavsky@aol.com" <beliavsky@127.0.0.1:7501> wrote:
>
>
>
>
> Hee. I forgot about that. And I don't think it's relevant for my
> question.
>
> Lapack/Scalapack is 500,000 lines of tremendously useful legacy code,
> and it's not going to go away anytime soon. At the moment my only
> ambition is to make that code as cleanly usable from other languages as
> possible. Using automatically translated versions is at the least a
> support nightmare.
>
> Thanks for your response,
>

The reasons that folks translate are more at the level of urban legends
and about as well informed as your wanting to retire Fortran.

The urban legends may relate to the problems of doing mixed languages.
You might require run time libraries for both but only have one
compiler. There is surely some amount of pure bigotry which may
be passed along by PHBs (Pointy Haired Boss - see Dilbert for a full
explanation). Policy may dictate that full source must be available
even if it is unreadable and unmaintainable in practice. This allows
the source to be recompiled on new and yet unknown machines. The new
machine may actually be the same old machine but with a new version
of its system software.

Each of these is of the "enduring nonsense" variety as there is a small
amount of technical justification which might apply is some circumstance
that typically has no relevance to the situation in which it is now
applied.

> V.

Victor Eijkhout

2004-12-10, 3:58 pm

Dan Pop <Dan.Pop@cern.ch> wrote:

> The issue that springs to mind is that Fortran and C store matrices
> differently: Fortran does it column-wise (the first index varies fastest),
> C does it row-wise (the last index varies fastest).
>
> The fact that Fortran array indexing is 1-based (by default) and C array
> indexing is 0-based shouldn't cause any problems when using the compiled
> Fortran code from C.

These sort of issues make me wonder why clapack is considered so great:
it does not provide arrays in the C-natural order, and in some places
the 1-based indexing is exposed, such as in pivoting of the
factorization.

When we write a more C-native interface we'll try to cover these issues.
It will mean some minimal bit of glue code that will hardly impact
performance.

Thanks for your response.

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Kamaraju Kusumanchi

2004-12-10, 3:58 pm

glen herrmannsfeldt

2004-12-10, 3:58 pm

Dan Pop wrote:

(snip regarding lapack and clapack)

> The issue that springs to mind is that Fortran and C store matrices
> differently: Fortran does it column-wise (the first index varies fastest),
> C does it row-wise (the last index varies fastest).

> The fact that Fortran array indexing is 1-based (by default) and C array
> indexing is 0-based shouldn't cause any problems when using the compiled
> Fortran code from C.

Some people who write numerical algorithms in Fortran exchange
the subscripts from the mathematical description, and others don't.

With cache memory and virtual memory one should always be sure that
the loops correspond to the memory layout in the appropriate language.

-- glen

Bill Shortall

2004-12-10, 8:58 pm

"Victor Eijkhout" <see.sig@for.addy> wrote in message
news:1gokv1s.16zyrkl1ytczkbN%see.sig@for.addy...
> Tino <tino52@yahoo.com> wrote:
>
>
> The thought has occurred to many people. I mean, passing LDA into a C++
> routine is the programming equivalent of turning a crank handle on the
> space shuttle. Now, what to do about it....
>
> I've started looking at Sidl/Babel to get a more native looking
> interface in other languages than F77. Not sure to what extent that will
> satisfy you.
>
> V.
> --
> email: lastname at cs utk edu
> homepage: www cs utk edu tilde lastname

Hi Victor,

For the last few years I have been working on a
set of C++ classes to do linear algebra. I was trying to
acheive most of the functionality of Lapack/Blas but
in a very user friendly fashion. The classes center
arround a general purpose vector class along with a dense rectangular Matrix
class and then some more compilcated classes like a vector of vectors and
vector of matrices.
All of them use operator overloading so you can write
A = B + C * D where A,B,C, D are matrices and all of the operators
work for both real and complex matrices. Akll of the basic Linear algebra
functions are supported
ie solution of equations, decompositions, SVD eigenvectors etc. etc.
The entire library is less than 5000 lines of code
and compiles into a static library of ~~ 3 megabytes
You don't have to compile it but you will need either a
Miccrosoft VC6 or Linux GCC compiler to use it.
I have demo version available and if anyone in the group wants to play
with it and has a compatible compiler send me a line. Once you've used it
you will never go back to LAPACK
pecos@cminet.net
regards...Bill

Victor Eijkhout

2004-12-10, 8:58 pm

Bill Shortall <pecos@cminet.net> wrote:

> For the last few years I have been working on a
> set of C++ classes to do linear algebra. I was trying to
> acheive most of the functionality of Lapack/Blas but
> in a very user friendly fashion. The classes center
> arround a general purpose vector class along with a dense rectangular Matrix
> class and then some more compilcated classes like a vector of vectors and
> vector of matrices.

Bill,

please mail me a pointer to your software.

However, without wanting to be disrespectful of your software, I wonder
if this is such a wise approach. Lapack is a high quality code base that
supports several data formats: dense, symmetric, banded, tridiagonal.
There is use for that, so before you can state

> Once you've used it
> you will never go back to LAPACK

you have to duplicate a lot of that functionality.

Then, Lapack has (especially in the eigenvalue part) state of the art
algorithms. The people coming up with the best algorithms are
implementing them straight in Lapack. What eigenvalue algorithms do you
use? Condition estimation, iterative refinement, et cetera.

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Michael Hosea

2004-12-10, 8:58 pm

I can't say whether my experience is rare or not, but there has not always
been a FORTRAN compiler available for the target platforms that I was
working on, namely the embedded ones. The C versions based on F2C weren't
very inviting, either. I can run F2C myself, and I have. I've even
hand-translated LAPACK code to efficient C. Even when a FORTRAN compiler is
available, it can be an expensive option when its only purpose is to compile
a FORTRAN library or two on a very occasional basis, and we might be talking
about more than one FORTRAN compiler if the C source were cross-compiled on
multiple platforms. In some cases it makes perfect economic sense to build
and use a FORTRAN library, and in others it doesn't seem as appealing.
Unfortunately for me, I've usually found myself in positions where it seemed
quite unappealing for one reason or another. On one system there was even a
custom C compiler that used a different floating point model. To my
knowledge, no compatible FORTRAN compiler ever even existed for that
platform!

What I have often wished for was a C implementation that represented the
same quality and attention to detail as the FORTRAN version.

I have no idea what you could do about any of that as long as FORTRAN is the
highest priority language to support, the "mother language", if you will.
Obviously FORTRAN is still heavily used in some niches. I'm not questioning
this, nor passing any kind of value judgement on it, but I would definitely
look very hard at the question of what the priority *should* be, basing this
on what you are trying to accomplish with a new version of LAPACK. You
might get one answer if, for example, you were primarily interested in
supercomputing or are only planning a relatively minor incremental update,
and you might reach a different conclusion if you wanted use of the library
to be easy and convenient to incorporate for the widest possible audience.
--
Mike

"Victor Eijkhout" <see.sig@for.addy> wrote in message
news:1gokq2z.1qtze21fbcwaoN%see.sig@for.addy...
> The authors of lapack/scalapack are starting to work on a new release of
> these packages. One of the things we want to address is the sillyness
> that C users, instead of linking to the binaries of the Fortran version,
> use an automatically translated C version.
>
> Therefore we'd like to know precisely what the reasons for this are.
>
> Naively one would say, compile the Fortran version, append an underscore
> (when appropriate) to routine names, and just link those libraries to
> your object files.
>
> What are the gotchas, and what are the real stumbling blocks here?
> Please spell it out in as much detail as you can muster. This situation
> really needs to improve.
>
> V.
> --
> email: lastname at cs utk edu
> homepage: www cs utk edu tilde lastname

Bill Shortall

2004-12-10, 8:58 pm

"Victor Eijkhout" <see.sig@for.addy> wrote in message
news:1gol29n.2hz9tynnu0bN%see.sig@for.addy...
> Bill Shortall <pecos@cminet.net> wrote:
>
Matrix[color=darkred]
and[color=darkred]
>
> Bill,
>
> please mail me a pointer to your software.
>
> However, without wanting to be disrespectful of your software, I wonder
> if this is such a wise approach. Lapack is a high quality code base that
> supports several data formats: dense, symmetric, banded, tridiagonal.
> There is use for that, so before you can state
>
>
> you have to duplicate a lot of that functionality.
>
> Then, Lapack has (especially in the eigenvalue part) state of the art
> algorithms. The people coming up with the best algorithms are
> implementing them straight in Lapack. What eigenvalue algorithms do you
> use? Condition estimation, iterative refinement, et cetera.
>
> V.
> --
> email: lastname at cs utk edu
> homepage: www cs utk edu tilde lastname

Hi Victor,
I need a e-mail address to send you the algebra package
its not on the internet
do you want the microsoft VC6 or the the Linux GCC 2.7
version. The package is about 750 kilobytes. I wonder
how big Lapack/Blas is ?
Sure the Lapack people might have some more sophisticated algorithms but
only a few problems need them. The advantages of working in C++ and having
all its tremendous libraries available to you are overwhelming.
The eigenvalue methods I use are the old standbys
reduction to either tri or bi diagonal then QR iteration. For
real general matrices I use reduction to hessian then QR iteration for the
eigenvalues and inverse iteration for the eigenvectors.

Victor Eijkhout

2004-12-10, 8:58 pm

Michael Hosea <mhosea@mathworks.com> wrote:

Very interesting response. Thanks.

I'm not sure that what we're doing right now is going to help you any.
However, Lapack has traditionally had a very clean coding style with
rigidly enforced style, so I imagine that translating automatically will
always stay a possibility.

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Mark McIntyre

2004-12-10, 8:58 pm

On Fri, 10 Dec 2004 16:59:32 -0500, in comp.lang.c , see.sig@for.addy
(Victor Eijkhout) wrote:

>Michael Hosea <mhosea@mathworks.com> wrote:

(stuff)

can you please stop crossposting to comp.lang.c as this is offtopic here.
thanks

--
Mark McIntyre
CLC FAQ <http://www.eskimo.com/~scs/C-faq/top.html>
CLC readme: <http://www.ungerhu.com/jxh/clc.welcome.txt>

Jentje Goslinga

2004-12-10, 8:58 pm

Victor Eijkhout wrote:
> The authors of lapack/scalapack are starting to work on a new release of
> these packages. One of the things we want to address is the sillyness
> that C users, instead of linking to the binaries of the Fortran version,
> use an automatically translated C version.
>
> Therefore we'd like to know precisely what the reasons for this are.

I am developing a suite of C functions consisting of the
Symmetric Eigenvalue problem, the SVD and the Block Reflector.
I have used Lapack as a reference to write a Divide and Conquer
solver for the Symmetric Eigenvalue problem and am now doing
the D&C SVD. I am quite well informed.

Many people nowadays perform numerical computations in C
or C++ which are just as suitable as Fortran.

There are many engineering applications with substantial
number crunching which are written in C.

Most individuals and many smaller Engineering companies do
not want to be perpetually burdened with having to purchase
a license for a Fortran compiler for their workstations and
would rather spend a few man days converting the Lapack
modules which their application requires to C.
Have you ever seen those bills for a Fortran compiler for an
AIX machine? Have you ever had to fight SUN Fortran?

Those mentioned would love to have the main factorizations in
C code rather than being burdened with having to purchase a
license for a Fortran compiler if one is at all available,
and having to cross language link yet another huge library
of which only a fraction is used.

I estimate there must be a fair contingent of people who
have extracted the CLapack modules required by their
application, found out what auxiliary modules are required
and absorbed these into a C application.

The problem with Lapack is that you can't read the code,
that it is neither commented nor documented with a Technical
Manual and that the algorithms referred to are detailed in
obscure internal reports at certain universities.
This is fine if there is reaaaallly no maintenance required
ever but this is never the case.
And Fortran to C conversion make the code even worse.
Although it is certainly possible to write code of almost the
same quality as good C in Fortran, the Lapack code uses 60's
style programming:
No dynamic memory allocation but the eternal WORK and IWORK
arrays which are partitioned up in intricate ways.
Does Lapack really need to cater to prehistorical compilers
which do not support dynamic memory allocation?
The functions generally are too large with the familiar gobs
of variables at the start and no attempt at resticting scope.
There are no Modules containing related functions but every
function is separate.
The single precision versions are useless in my opinion but
others may disagree, at the very least relegate them to
recognizeable modules so users who do not want them can leave
them out.
There is more: insightful remarks on how to re-use arrays to
conserve memory, fear of copying arrays resulting in a milliped
matrix multiply functionality and so on.
"Lame comments" and huge headers abound but the intricacies of
the algorithms are not mentioned.

If you want to sink more time into the Fortran Lapack, you
might want to look into the temporary storage issue, allocating
the memory instead so the user is not burdened with having to
calculate the "work storage" required.
Next get rid of the single precision versions or at least
modularize them and generally aggregate the code into a few
dozen modules of related functions.
For the QR type methods, store the Givens transformations and
allow the user to pass a function to select Eigenvectors or
Singular Vectors for computation rather than having to estimate
their index for computation beforehand which really sucks.
Just some pointers.

Lapack's use of a BLAS is great but Lapack is tailored to the
traditional BLAS rather than the other way around and the
functionality of these BLAS functions is not optimal.
Simple Example: BLAS provides DNORM2 and DSCALE functions but
90% plus of the use of these is to Unitize a vector which should
be a single function. Currently this requires two function calls
between which the all important norm information is truncated
from 80 bits accuracy to 64 on a PC for example. Likewise with
the computation of the Householder vector. Is this a big deal?
Well, there are quite a few of these instances, and in the end
it all adds up. The point is that you might want to look beyond
the traditional BLAS which is just historical accident.

Lapack is quite accurate and fast although the accuracy of some
of the order two factorizations and some of the orthogonal
transformations can still be improved here and there.

> Naively one would say, compile the Fortran version, append an underscore
> (when appropriate) to routine names, and just link those libraries to
> your object files.
>
> What are the gotchas, and what are the real stumbling blocks here?
> Please spell it out in as much detail as you can muster. This situation
> really needs to improve.

I have spelled it out; this is not intended as an attack but
as constructive criticism.
Lapack is free and has been useful to me as a reference and
for accuracy tsting my own code.
I realize that a lot of people have spent considerable effort
in its design and I am grateful for that. And it works.

Thanks,

> V.

Jentje Goslinga

my home site is:
http://www.numerical-algorithms.com/

Greg Lindahl

2004-12-10, 8:58 pm

In article <41BA2BB2.6060003@telus.net>,
Jentje Goslinga <goslinga@telus.net> wrote:

> Currently this requires two function calls
> between which the all important norm information is truncated
> from 80 bits accuracy to 64 on a PC for example.

Wow. I can't believe anyone would post this to sci.math.num-analysis!
I thought that non-explicit use of extended precision variables was
universally considered to be a dumb idea. If you want explicit use
of extended precision with LAPACK, you need extra functions.

BTW, modern compilers on PCs generally default to using SSE2
instructions on modern PC cpus, hence no more extended precision
unless you ask for it. This became true soon after the Pentium 4 was
released, which was a few years ago.

-- greg

Michael Hosea

2004-12-10, 8:58 pm

This is what I love about usenet. You can respond to one question every 4
years or so, and you immediately get asked to "stop". :-)
--
Mike

"Mark McIntyre" <markmcintyre@spamcop.net> wrote in message
news:bl9kr0dm6la8r0uusjc27fa450ucq3n72u@
4ax.com...
> On Fri, 10 Dec 2004 16:59:32 -0500, in comp.lang.c , see.sig@for.addy
> (Victor Eijkhout) wrote:
>
>
> (stuff)
>
> can you please stop crossposting to comp.lang.c as this is offtopic here.
> thanks
>
>
> --
> Mark McIntyre
> CLC FAQ <http://www.eskimo.com/~scs/C-faq/top.html>
> CLC readme: <http://www.ungerhu.com/jxh/clc.welcome.txt>

beliavsky@aol.com

2004-12-11, 3:56 am

>Many people nowadays perform numerical computations in C
>or C++ which are just as suitable as Fortran.

Fortran 95 still has many advantages over C and C++:
(1) true multidimensional arrays, with a wide range of intrinsic
functions like SUM and MAXVAL and the ability to use array sections
like a(1:10,:,2:10:2).
(2) the ability to pass one- and multi-dimensional arrays to functions
in a SIMPLE manner and to have functions returning arrays.
(3) arrays that can start at 0, 1, or any other user-specifed integer
(4) intrinsic functions like SIN that are generic to single, double,
and often quadruple precision, for both scalar and array arguments
(5) user-defined PURE and ELEMENTAL functions, to facilitate
parallelization and reduce the need for separate functions handling
scalar and array arguments
(6) better handling of complex numbers
(7) a form of DO loop that prevents the user from inadvertently
changing the loop variable. WHILE loops and and a DO-ENDDO loop with an
EXIT can be used when more flexibility is needed.
(8) the ability to EXIT from nested loops
superior facilities for reading and writing formatted tables of numbers
(9) a cleaner, less error-prone syntax. It is natural for the
end-of-line to terminate a statement by default, and to allow
continuation lines when needed. A misplaced semi-colon won't
drastically change the meaning of a Fortran code.

>There are many engineering applications with substantial
>number crunching which are written in C.

In my company substantial number crunching is done in Excel, because
some people refuse to learn anything else. Many C programmers seem to
have similar blinders.

>Most individuals and many smaller Engineering companies do
>not want to be perpetually burdened with having to purchase
>a license for a Fortran compiler for their workstations and
>would rather spend a few man days converting the Lapack
>modules which their application requires to C.
>Have you ever seen those bills for a Fortran compiler for an
>AIX machine? Have you ever had to fight SUN Fortran?

Much scientific computation that used to be done on commercial UNIX
platforms is now being done on Linux, where there are many commercial
Fortran 95 compilers. The Lahey/Fujitsu Linux compiler costs as little
as $250, and the Intel Fortran compiler is free for non-commercial use.
G95 at http://www.g95.org is 99.9% of a free, open-source Fortran 95
compiler for Linux, Windows, Mac OS X, and other platforms. Gfortran is
getting there.

>The problem with Lapack is that you can't read the code,
>that it is neither commented nor documented with a Technical
>Manual and that the algorithms referred to are detailed in
>obscure internal reports at certain universities.

Lapack documentation is at http://www.netlib.org/lapack/lug/index.html
, also as a printed book.

>the Lapack code uses 60's style programming:
>No dynamic memory allocation but the eternal WORK and IWORK
>arrays which are partitioned up in intricate ways.
>Does Lapack really need to cater to prehistorical compilers
>which do not support dynamic memory allocation?

Fortran 90/95 has dynamic memory allocation.

>The functions generally are too large with the familiar gobs
>of variables at the start and no attempt at resticting scope.

Fortran 90/95 functions and subroutines can have OPTIONAL arguments, so
that the user may only need to specify a few to get the default
behavior. The LAPACK95 interface to the LAPACK library already exists,
and using it, a subroutine call such as

call SGELSD(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK,LWORK,
IWORK, INFO)

can be replaced with just

CALL LA_GELSD (A, B)

>There are no Modules containing related functions but every function
is separate.

Fortran 90/95 has something called a "MODULE" that does exactly this.
To get compile-time checking of procedure interfaces, the
"pre-historic" #include statements of C/C++ are not needed.

>The single precision versions are useless in my opinion but
>others may disagree, at the very least relegate them to
>recognizeable modules so users who do not want them can leave
>them out.

The single precision codes of LAPACK start with the letter S, as
described at http://www.netlib.org/lapack/individualroutines.html , and
are thus easy to recognize.

Ideally, each of the BLAS calls could be replaced by calls to Fortran
90/95 intrinsic functions or array operations, improving readability.

You have indirectly made a good case for writing the next version of
LAPACK in Fortran 95.

In general, one should try to gain a current understanding of a
subject, in this case Fortran, before making recommendations about it
to others. There are Fortran 90/95 tutorials at
http://www.dmoz.org/Computers/Progr...tran_90_and_95/
..

James Van Buskirk

2004-12-11, 3:56 am

<beliavsky@aol.com> wrote in message
news:1102727931.346615.98970@c13g2000cwb.googlegroups.com...

> The single precision codes of LAPACK start with the letter S, as
> described at http://www.netlib.org/lapack/individualroutines.html , and
> are thus easy to recognize.

What I would prefer regarding the simgle precision versions would
be a single unified version, with numeric inquiry functions, where
necessary, accounting for the differences. The kind type parameter
in the functions would be inherited from the host module. Another
module could rename the procedures to their single- and double-
precision names and create also generic names. The utility of the
resulting code would be that of C++ templates in that all the code
would be maintained in one place and could be used generically. The
only comparative awkwardness of this approach is the intermediate
module(s) that do the renaming, but I don't think they would be
changed a whole lot in the maintainance phase. I have previously
posted several examples of this technique, but it doesn't seem to
have caught on.

> You have indirectly made a good case for writing the next version of
> LAPACK in Fortran 95.

My reading of the previous post was that the intent was more direct.
But maybe I'm just reading others' posts more optimistically.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

Allin Cottrell

2004-12-11, 3:56 am

Dan Pop wrote:
> In <1gokq2z.1qtze21fbcwaoN%see.sig@for.addy> see.sig@for.addy (Victor Eijkhout) writes:
>
>
>
>
> The issue that springs to mind is that Fortran and C store matrices
> differently: Fortran does it column-wise (the first index varies fastest),
> C does it row-wise (the last index varies fastest).

I don't think that is the issue, since the clapack translation (as I
understand it) is not an idiomatic translation into C (actually, that
would be much too confusing, if you think about it). It's a perfectly
literal translation, which preserves Fortan-style array indexing by
brutal mangling of C notation.

Allin Cottrell
Wake Forest University

Jentje Goslinga

2004-12-11, 3:56 am

beliavsky@aol.com wrote:
>
>
> Fortran 95 still has many advantages over C and C++:
> (1) true multidimensional arrays, with a wide range of intrinsic
> functions like SUM and MAXVAL and the ability to use array sections
> like a(1:10,:,2:10:2).

That is great, now tell me why the part of Lapack, at least the
functions for Eigenvalues and SVD I have looked at does not use
any of these.

> (2) the ability to pass one- and multi-dimensional arrays to functions
> in a SIMPLE manner and to have functions returning arrays.
> (3) arrays that can start at 0, 1, or any other user-specifed integer
> (4) intrinsic functions like SIN that are generic to single, double,
> and often quadruple precision, for both scalar and array arguments
> (5) user-defined PURE and ELEMENTAL functions, to facilitate
> parallelization and reduce the need for separate functions handling
> scalar and array arguments
> (6) better handling of complex numbers

That is great but I rarely need those and can write the required
operations in an hour.

> (7) a form of DO loop that prevents the user from inadvertently
> changing the loop variable. WHILE loops and and a DO-ENDDO loop with an
> EXIT can be used when more flexibility is needed.

Come on, I am talking to programmers.

> (8) the ability to EXIT from nested loops
> superior facilities for reading and writing formatted tables of numbers
> (9) a cleaner, less error-prone syntax. It is natural for the
> end-of-line to terminate a statement by default, and to allow
> continuation lines when needed. A misplaced semi-colon won't
> drastically change the meaning of a Fortran code.

Sigh.

>
>
> In my company substantial number crunching is done in Excel, because
> some people refuse to learn anything else. Many C programmers seem to
> have similar blinders.

I never called myself exclusively a C programmer and may have
written a lot more Fortran code than you have.

>
>
> Much scientific computation that used to be done on commercial UNIX
> platforms is now being done on Linux, where there are many commercial
> Fortran 95 compilers. The Lahey/Fujitsu Linux compiler costs as little
> as $250, and the Intel Fortran compiler is free for non-commercial use.
> G95 at http://www.g95.org is 99.9% of a free, open-source Fortran 95
> compiler for Linux, Windows, Mac OS X, and other platforms. Gfortran is
> getting there.

It all depends on the application. Nowadays "systems" usually
have a database component and sometimes internet connectivity
and much more and it is natural for them to be composed of
modules in various different languages. I have seen many
engineering applications in the C and C++ languages.
I have written a large system with a Windows user interface,
database and internet connectivity in C++.
The beauty of C++ is that the notion that the code will have
to connect to modules written in other languages has been an
intrinsic part of the design of the language from the beginning.

>
>
> Lapack documentation is at http://www.netlib.org/lapack/lug/index.html
> , also as a printed book.

I haven't looked there but assume these are user manuals
for people like you. I am talking about design here.

>
>
> Fortran 90/95 has dynamic memory allocation.

Precisely, that is what my post implies.

>
>
> Fortran 90/95 functions and subroutines can have OPTIONAL arguments, so
> that the user may only need to specify a few to get the default
> behavior. The LAPACK95 interface to the LAPACK library already exists,
> and using it, a subroutine call such as
>
> call SGELSD(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK,LWORK,
> IWORK, INFO)
>
> can be replaced with just
>
> CALL LA_GELSD (A, B)
>
>
>
> is separate.
>
> Fortran 90/95 has something called a "MODULE" that does exactly this.
> To get compile-time checking of procedure interfaces, the
> "pre-historic" #include statements of C/C++ are not needed.
>
>
>
>
> The single precision codes of LAPACK start with the letter S, as
> described at http://www.netlib.org/lapack/individualroutines.html , and
> are thus easy to recognize.

Yes, I know that. There is rarely an excuse for duplicating code
such as for different precisions for example, this is really a
design problem. Beware... before you make that copy.

> Ideally, each of the BLAS calls could be replaced by calls to Fortran
> 90/95 intrinsic functions or array operations, improving readability.
>
> You have indirectly made a good case for writing the next version of
> LAPACK in Fortran 95.

I am making lots of useful and constructive comments.

> In general, one should try to gain a current understanding of a
> subject, in this case Fortran, before making recommendations about it
> to others.

Exactly, that's why I sometimes take the liberty to post.
I have programmed in Fortran since 1980, in C since 1985 and
in C++ and Assembler since 1995.

There are Fortran 90/95 tutorials at
>
http://www.dmoz.org/Computers/Progr...tran_90_and_95/

Enough!

The purpose of my post was not to discredit Fortran and
I never did, this is not another language war.
My article makes some genuinely useful and well researched
suggestions for improvement which are based on a lot of
thinking. You may think otherwise.

Jentje Goslinga

beliavsky@aol.com

2004-12-11, 3:56 am

Jentje Goslinga <goslinga@telus.net> wrote:
>
>
>beliavsky@aol.com wrote:
>
>That is great, now tell me why the part of Lapack, at least the
>functions for Eigenvalues and SVD I have looked at does not use
>any of these.

Because it is written in Fortran 77, which did not have these features. A
new LAPACK should be in Fortran 95 IMO.

----== Posted via Newsfeeds.Com - Unlimited-Uncensored-Secure Usenet News==----
http://www.newsfeeds.com The #1 Newsgroup Service in the World! >100,000 Newsgroups
---= 19 East/West-Coast Specialized Servers - Total Privacy via Encryption =---

Rich Townsend

2004-12-11, 3:56 am

Jentje Goslinga wrote:

> It all depends on the application. Nowadays "systems" usually
> have a database component and sometimes internet connectivity
> and much more and it is natural for them to be composed of
> modules in various different languages. I have seen many
> engineering applications in the C and C++ languages.
> I have written a large system with a Windows user interface,
> database and internet connectivity in C++.
> The beauty of C++ is that the notion that the code will have
> to connect to modules written in other languages has been an
> intrinsic part of the design of the language from the beginning.

I can't find the section in Bjarne Stroustrup's book where he discusses
interoperability with other languages. Would you care to enlighten me?

cheers,

Rich

--
Dr Richard H D Townsend
Bartol Research Institute
University of Delaware

[ Delete VOID for valid email address ]

Jentje Goslinga

2004-12-11, 3:56 am

Rich Townsend wrote:
> Jentje Goslinga wrote:
>
>
>
> I can't find the section in Bjarne Stroustrup's book where he discusses
> interoperability with other languages. Would you care to enlighten me?

Yes, the "extern C" syntax is discused in Chapter 9
Section 9.2.4 Linkage to Non-C++ code,
The C++ programming Language, third edition page 205.

There is also a useful discussion about the design of large
systems in Chapter 23, Development and Design, worth reading.

This a great book.

> cheers,

Thanks for saying something positive for a change.

Jentje

> Rich

James Giles

2004-12-11, 3:56 am

Jentje Goslinga wrote:
> Rich Townsend wrote:
....[color=darkred]
>
> Yes, the "extern C" syntax is discused in Chapter 9
> Section 9.2.4 Linkage to Non-C++ code,
> The C++ programming Language, third edition page 205.

That's interop with C (which is almost redundant, C++ originally
being an extension of C, it's hardly impressive that it can interop
with C). I still don't see the part where it's originally intended
to interop with, say, Ada, Pascal, Lisp, Haskell, and yes: Fortran.

--
J. Giles

"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare

Gerry Thomas

2004-12-11, 3:56 am

"Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
news:cpds8l$kop$1@scrotar.nss.udel.edu...
>
> I can't find the section in Bjarne Stroustrup's book where he discusses
> interoperability with other languages. Would you care to enlighten me?
>

Duh!, dummy Andy Capp from U of Delaware, with C++, who needs Fortran or
any other language?. There are several Fortran compilers (sorry, you have
to buy them, a known problem for cheapskates like you) that work quite
harmoniously with C/C++ but you're not ready for them, you being way too
dimwitted.
--
You're Welcome,
Gerry T.
______
"Ah, Klinger, my constant reminder that Darwin was right!." -- Maj Charles
Winchester III, the 4077th M*A*S*H

Rich Townsend

2004-12-11, 3:56 am

Jentje Goslinga wrote:
>
>
> Rich Townsend wrote:
>
>
>
> Yes, the "extern C" syntax is discused in Chapter 9
> Section 9.2.4 Linkage to Non-C++ code,
> The C++ programming Language, third edition page 205.

Page 524 (r.7.4.) in my second edition (I assume we are talking about
"The C++ Programming Language"). But in all seriousness, you can't claim
that 'extern' is in any way comparable with the sort of language
operability we see nowadays. The meaning of the string literal following
the 'extern' is implementation dependent, therefore this mechanism for
interoperability is wholly non-portable. Fortran has has similar
non-portable language interoperability *long* before C++ even existed --
there are few mainstream f77 compilers that aren't interoperable with C
on a platform-dependent basis.

Therefore, your claim about C++ being designed ab initio with
interoperability in mind is a little misleading; an equally-meaningless
claim could be made of Fortran, of Pascal, of Forth, of Basic etc etc.

>
> There is also a useful discussion about the design of large
> systems in Chapter 23, Development and Design, worth reading.
>
> This a great book.

Indeed; but I find the style a little dry. Metcalf & Reid is closer to
my preferred pedagogical style -- not as flippant as Numerical Recipes,
not as dense as Knuth.

>
>
>
> Thanks for saying something positive for a change.

Out of curiousity, me in particular, or just people in general?

cheers,

Rich

--
Dr Richard H D Townsend
Bartol Research Institute
University of Delaware

[ Delete VOID for valid email address ]

Victor Eijkhout

2004-12-11, 3:56 am

Jentje Goslinga <goslinga@telus.net> wrote:

> No dynamic memory allocation but the eternal WORK and IWORK
> arrays which are partitioned up in intricate ways.
> Does Lapack really need to cater to prehistorical compilers
> which do not support dynamic memory allocation?

Ok, it's clear where this use of temporaries comes from. However, it's
not clear that dynamic allocation is the best solution. If your routine
is called once, with a large dataset, then it could do its own
alloation. However, if the routine is called a large number of times
with the same small problem size, then you want to do the allocation
outside it.

We probably have to take a step back to the application level, and see
which routines can be part of a computationally intensive inner loop.

Your comments very much appreciated.

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Victor Eijkhout

2004-12-11, 3:56 am

"beliavsky@aol.com" <beliavsky@127.0.0.1:7501> wrote:

> Because it is written in Fortran 77, which did not have these features. A
> new LAPACK should be in Fortran 95 IMO.

I'm sorry, that's too much code to rewrite.

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Jentje Goslinga

2004-12-11, 3:56 am

James Giles wrote:
> Jentje Goslinga wrote:
>
> ...
>
>
>
>
> That's interop with C (which is almost redundant, C++ originally
> being an extension of C, it's hardly impressive that it can interop
> with C).

It is not so trivial with overloading and name decoration.

> I still don't see the part where it's originally intended
> to interop with, say, Ada, Pascal, Lisp, Haskell, and yes: Fortran.

You are right, that part doesnt exist; I had forgotten that
but what's the difference, isn't C almost the same as Fortran?

Interoperability between languages is a gray area where one
needs to consider the following:
1 Ability to suppress C++ name decoration, knowledge of other
conventions such as underscores appended or prepended to name
2 Knowledge of linker conventions such as upper case conversion
and maximum name length
3 Ability for the programmer to specify parameter passing by
address and by value which C and C++ have
4 Ability for the user to use parameter passing from left to
right or vice versa
5 Conventions on stack based parameter passing notably padding
and stack alignment
6 Decide on who resets the stack pointer, caller or callee.

I agree there may be differences between languages between
items 4-6 although they are usually not irreconcilable since
the operating systems supporting these languages must use
cross language linking too.
I believe none of these are even specified by the C or C++
standards either, although parameter passing from the right
to left is said to be necessary to accommodate variable
number of parameters.

But you are right.

Jentje

Jentje Goslinga

2004-12-11, 3:56 am

Rich Townsend wrote:
> Jentje Goslinga wrote:
>
>
>
> Page 524 (r.7.4.) in my second edition (I assume we are talking about
> "The C++ Programming Language"). But in all seriousness, you can't claim
> that 'extern' is in any way comparable with the sort of language
> operability we see nowadays. The meaning of the string literal following
> the 'extern' is implementation dependent, therefore this mechanism for
> interoperability is wholly non-portable. Fortran has has similar
> non-portable language interoperability *long* before C++ even existed --
> there are few mainstream f77 compilers that aren't interoperable with C
> on a platform-dependent basis.
>
> Therefore, your claim about C++ being designed ab initio with
> interoperability in mind is a little misleading; an equally-meaningless
> claim could be made of Fortran, of Pascal, of Forth, of Basic etc etc.
>

You have said it.

>
>
> Indeed; but I find the style a little dry. Metcalf & Reid is closer to
> my preferred pedagogical style -- not as flippant as Numerical Recipes,
> not as dense as Knuth.
>
>
>
> Out of curiousity, me in particular, or just people in general?

No, you, thanks for saying something positive. This newsgroup
stuff can be quite draining after a while especially if one's
expertise is continually questioned.
I would not post if I didnt think I had something to contribute
based on experience.

Jentje

> cheers,
>
> Rich
>

James Giles

2004-12-11, 3:56 am

Jentje Goslinga wrote:
> James Giles wrote:
....
>
> You are right, that part doesnt exist; I had forgotten that
> but what's the difference, isn't C almost the same as Fortran?

Well, not really. Historically is was not at all infrequent to find
that many languages could not interoperate easily with C, even
though they could easily interoperate with each-other. C often
tended to be diffferent, sometimes for no discerable reason: just
to be different I guess.

Also you mentioned "name decoration" or "name mangling" as
it's more often called. (It is, of course, the spawn of the devil.)
That's one of the cases in which many other languages could
interoperate, and no "mangling" was needed, while C
implementations, again often for no real reason, insisted
that non-C names *had* to be mangled. I never saw name
mangling (and I did site support of production linkers) before
C entered general service.

Note: I understand the reason for name mangling, but that doesn't
mean I have to like it. It's a lazy implementation choice that shifts
some of the burden elsewhere, including the users, when it should
be managed more cleanly by the compilers and linkers.

--
J. Giles

"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare

Jentje Goslinga

2004-12-11, 3:56 am

Victor Eijkhout wrote:
> Jentje Goslinga <goslinga@telus.net> wrote:
>
>
>
>
> Ok, it's clear where this use of temporaries comes from. However, it's
> not clear that dynamic allocation is the best solution. If your routine
> is called once, with a large dataset, then it could do its own
> alloation. However, if the routine is called a large number of times
> with the same small problem size, then you want to do the allocation
> outside it.

That is very true but in many cases it may be moved a level
below the user.

> We probably have to take a step back to the application level, and see
> which routines can be part of a computationally intensive inner loop.
>
> Your comments very much appreciated.

Thank you for saying that; I am trying to be constructive
and freely admit to having benefitted from the availability
of Lapack during the testing of my own code.

I also agree with you that it would be desirable to integrate
the single/double precision functions, if possible, preferably
prior to the 16 byte reals appearing on the scene.

> V.

Jentje

glen herrmannsfeldt

2004-12-11, 8:56 am

Greg Lindahl wrote:
(snip)

> Wow. I can't believe anyone would post this to sci.math.num-analysis!
> I thought that non-explicit use of extended precision variables was
> universally considered to be a dumb idea. If you want explicit use
> of extended precision with LAPACK, you need extra functions.

> BTW, modern compilers on PCs generally default to using SSE2
> instructions on modern PC cpus, hence no more extended precision
> unless you ask for it. This became true soon after the Pentium 4 was
> released, which was a few years ago.

As far as I know, about half the time the extra precision helps,
and about half the time it hurts. (It hurts when things
round the wrong way, or when quantities that should be equal
are not equal.)

-- glen

Salvatore

2004-12-11, 3:58 pm

Victor Eijkhout wrote:
> Ok, it's clear where this use of temporaries comes from. However,
it's
> not clear that dynamic allocation is the best solution. If your
routine
> is called once, with a large dataset, then it could do its own
> alloation. However, if the routine is called a large number of times
> with the same small problem size, then you want to do the allocation
> outside it.
>
> We probably have to take a step back to the application level, and
see
> which routines can be part of a computationally intensive inner loop.
>

If you do write a shell in Fortran 95 it would be easy to have an
optional WORK parameter that allows the user to do his/her own
allocation when needed, and the shell will allocate its own otherwise
before calling the inner code (be it legacy F77 or newly rewritten F95
or whatever).

Cheers
Salvatore

beliavsky@aol.com

2004-12-11, 8:56 pm

> Jentje Goslinga <gosli...@telus.net> wrote:

<snip>

>I also agree with you that it would be desirable to integrate
>the single/double precision functions, if possible, preferably
>prior to the 16 byte reals appearing on the scene.

16-byte reals have already appeared on the scene. Several Fortran 95
compilers will compile the following code, which declares 4, 8, and
16-byte real variables. (It would be better to use the
SELECTED_REAL_KIND intrinsic function.)

real (kind=4) :: xs
real (kind=8) :: xd
real (kind=16) :: xq

REAL and DOUBLE PRECISION variables are usually stored in 4 and 8
bytes, but I believe the standard does not require this. The Numerical
Recipes in Fortran 90 library does not declare variables as REAL or
DOUBLE PRECISION but uses the KIND specifier instead. A new LAPACK
probably ought to do the same, which is another reason to write it in
Fortran 95.

glen herrmannsfeldt

2004-12-11, 8:56 pm

beliavsky@aol.com wrote:
(snip)

> 16-byte reals have already appeared on the scene. Several Fortran 95
> compilers will compile the following code, which declares 4, 8, and
> 16-byte real variables. (It would be better to use the
> SELECTED_REAL_KIND intrinsic function.)

At least since the 360/85 in about 1968. Also, CDC machines
with 60 bit words had 120 bit double precision.

> real (kind=4) :: xs
> real (kind=8) :: xd
> real (kind=16) :: xq

> REAL and DOUBLE PRECISION variables are usually stored in 4 and 8
> bytes, but I believe the standard does not require this. The Numerical
> Recipes in Fortran 90 library does not declare variables as REAL or
> DOUBLE PRECISION but uses the KIND specifier instead. A new LAPACK
> probably ought to do the same, which is another reason to write it in
> Fortran 95.

Well, the standard tends to require that default REAL have
the same size as default INTEGER, and that DOUBLE PRECISION
have twice that size. (I don't believe it requires all the
bits to actually be used, but the storage used must agree.)

-- glen

Victor Eijkhout

2004-12-11, 8:56 pm

Salvatore <sfilippone@uniroma2.it> wrote:

> If you do write a shell in Fortran 95 it would be easy to have an
> optional WORK parameter that allows the user to do his/her own
> allocation when needed, and the shell will allocate its own otherwise
> before calling the inner code (be it legacy F77 or newly rewritten F95
> or whatever).

That's a

idea. Thanks for the suggestion.

(Hope your trip home was uneventful?)

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Robert Corbett

2004-12-12, 3:58 am

In article <41BA2BB2.6060003@telus.net>,
Jentje Goslinga <goslinga@telus.net> wrote:
>
Have you ever had to fight SUN Fortran?

I am very interested in learning what problems you
are having with Sun Fortran. Please let me know what
they are. Please also let me know which release of
Sun Fortran you are using.

Yours truly,
Robert Corbett

hchu@thermawave.com

2004-12-12, 3:58 am

>
> Ok, it's clear where this use of temporaries comes from. However,
it's
> not clear that dynamic allocation is the best solution. If your
routine
> is called once, with a large dataset, then it could do its own
> alloation. However, if the routine is called a large number of times
> with the same small problem size, then you want to do the allocation
> outside it.
Nowadays, I typically work with array dimensions between 10 and 30
I have done a lot of tests and noticed that there is no running time
difference whether you repeatly use the same work arrays or allocate
them on demand.

In one of the least squares routines, the work array dimension is
very difficult to figure out. As other people suggested it would be
very nice to have some wrapper routines available.

We use CVF and Intel Fortran on both Windows and Linux, one of the
most annoying thing related to the lapack routines is the string
arguments. Because we use both C++ and Fortran and there is no
standard for string implementations, I have to write wrapper routines
to deal with that, like passing integers.

Another complaint I have is that a few routines are not thread safe
because some variables have to be computed and saved. It took me quite
some time to search for the SAVE statements and eventually wrote a
routine to call these routines before threads are created. It would
be nice for LAPACK to provide such routine.

Victor Eijkhout

2004-12-12, 3:58 am

<hchu@thermawave.com> wrote:

Thanks for your response.

> It took me quite
> some time to search for the SAVE statements and eventually wrote a
> routine to call these routines before threads are created. It would
> be nice for LAPACK to provide such routine.

Yeah, SAVE statements are on the "must get rid of" list.

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

glen herrmannsfeldt

2004-12-12, 3:58 am

Victor Eijkhout wrote:

> <hchu@thermawave.com> wrote:

[color=darkred]
> Yeah, SAVE statements are on the "must get rid of" list.

As far as I know, compilers are allowed to SAVE data,
even without a SAVE statement. Traditionally, Fortran
used static storage, )though with some complications due to
overlays) and wouldn't be appropriate for multithreaded code.

But yes, SAVE should only be used for data that really needs
to be static, a favorite example being the seed for a
random number generator.

-- glen

Jentje Goslinga

2004-12-12, 3:58 am

Robert Corbett wrote:
> In article <41BA2BB2.6060003@telus.net>,
> Jentje Goslinga <goslinga@telus.net> wrote:
>
> Have you ever had to fight SUN Fortran?
>
> I am very interested in learning what problems you
> are having with Sun Fortran. Please let me know what
> they are. Please also let me know which release of
> Sun Fortran you are using.
>
> Yours truly,
> Robert Corbett

Well, if you really ask for it: I have used Sun C and
Fortran from 1985-1990.
The transcendentals were only ten decimals accurate and
reported bugs were not fixed but maybe that was not in
Fortran. The whole experience with the hardware and the
software was traumatic throughout.

In 2000 I had the misfortune to have to work Sun machines
again, this time in C.
The simulation model I was using gave quite different
results on different machines with different versions
of the libraries. I simply don't trust Sun floating point.

Sun likes to call their boxes "workstations" but I
consider an Intel PC a much superior workstation for
floating point computations, and a lot cheaper.

Jentje Goslinga

Ron Shepard

2004-12-12, 3:58 pm

In article <slrncroaij.tik.romanNOSPAM@theta1.cft.edu.pl>,
Roman Werpachowski <"r o m a nNOSPAM"@theta1.cft.edu.pl> wrote:

> The only sensible way would be to have two functions: one which does the
> allocation, another which does not. The problem is, that sometimes one
> can flexibly choose a work space area, like in (D,S)SYEV routine. Small
> workspace conserves memory but large workspace gives better convergence.
> Which options should the self-allocating version of (D,S)SYEV choose?

Modern Fortran allows optional arguments, so the "sensible" way in
that language would be to specify the workspace with an optional
array. If specified during the call, then it is used, otherwise the
subroutine allocates its own workspace. Of course, this puts
additional burdens on other languages in cross-language projects to
support the fortran calling conventions. BTW, fortran also allows
the call to use a generic name (e.g. call syev() in the above
example), which is resolved to the specific precision based on the
data types of its arguments, in order to eliminate the (D,S)
portability problems when moving code between machines.

$.02 -Ron Shepard

Jentje Goslinga

2004-12-12, 8:57 pm

Bill Shortall wrote:
> "Victor Eijkhout" <see.sig@for.addy> wrote in message
> news:1gokv1s.16zyrkl1ytczkbN%see.sig@for.addy...
>
>
>
> Hi Victor,
>
> For the last few years I have been working on a
> set of C++ classes to do linear algebra. I was trying to
> acheive most of the functionality of Lapack/Blas but
> in a very user friendly fashion. The classes center
> arround a general purpose vector class along with a dense rectangular Matrix
> class and then some more compilcated classes like a vector of vectors and
> vector of matrices.
> All of them use operator overloading so you can write
> A = B + C * D where A,B,C, D are matrices and all of the operators
> work for both real and complex matrices.

This is not a very common operation which I have only
encountered when working with block matrices.
C++ is a powerful language but the simple algebra of
Matrices, vectors and Scalars is sufficiently complex
that using C++ shorthand notation is not always possible.
Have you looked at OON, the Object Oriented Numerics
Group?

Akll of the basic Linear algebra
> functions are supported
> ie solution of equations, decompositions, SVD eigenvectors etc. etc.
> The entire library is less than 5000 lines of code

You must be kidding, by the time I am done with Schur, SVD
and Block Reflector, I will have 25-30k lines of C and a few
thousand lines of Intel Assembler. And that is without Linear
Equations or Complex Data type.
And it is including testing modules for everything.

> and compiles into a static library of ~~ 3 megabytes

Isn't that a bit large? My executable which contains all of
the said factorizations and comprehensive testing software
weighs in at 76 kilobytes. You generate 40 times the object
code with one fifth the code, a factor 200!

> You don't have to compile it but you will need either a
> Miccrosoft VC6 or Linux GCC compiler to use it.
> I have demo version available and if anyone in the group wants to play
> with it and has a compatible compiler send me a line. Once you've used it
> you will never go back to LAPACK
> pecos@cminet.net
> regards...Bill

Have you tested those functions by constructing pathological
matrices from their eigenvalues (for symmetric matrices)
or from their singular values (in the case of the SVD) and
likewise for Linear Equations and evaluated the accuracy?

Jentje Goslinga

James Giles

2004-12-12, 8:57 pm

Well, I agree and disagree to some extent with both sides in this
thread so far, so I've mostly been watching. But, there are so
few absolutes in computing that I kind of shudder when I see
invalid ones promoted.

Jentje Goslinga wrote:
....
> Functions should be self-contained and allocate local storage.
> If they really perform so little actual work that this is a
> concern, then maybe the program is not structured optimally.

Maybe, maybe not. Define "optimal structure". Lots of things
are structured better when procedures do little work. Whatever
makes the code more legible and easier to maintain, is usually
quite superior. Contrary to the state of things in the '50s, human
time is usually now more valuable than the computer's time.
If allocating workspace outside of the procedures retains
the advantages of legibility while improving computer efficiency,
it's a valid technique. (Nor can I think of any way to prohibit
the technique with a language design change. So this is independent
of the language you use.)

> Variable arguments are another dubious language feature which
> made its way into the C language because of the requirement of
> the formatted print functions, but it causes problems.

C's feature is abominably defined. But, as long as you permit
generic procedures, it's hard to prohibit Fortran's form of optional
arguments. Fortran just has a notation that makes them more
convenient to write than as generics. But, from the point of
view of the user of the procedure, F(X, Other) vs. F(X) is the
same whether F is really two procedures that are generic (have
the same name, but different type signatures) or a single procedure
with an optional argument. (So, once again, you could actually
use this trechnique independent of what language you use.)

--
J. Giles

"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare

Mario S. Mommer

2004-12-12, 8:57 pm

Jentje Goslinga <goslinga@telus.net> writes:
> Variable arguments are another dubious language feature which
> made its way into the C language because of the requirement of
> the formatted print functions, but it causes problems.

Mostly because the way C "supports" such a thing is so ridiculously
baroque.

> Esoteric language "features" like variable arguments and "multiple
> returns" should be avoided if at all possible.

Variable arguments, keyword arguments, and multple value returns are
extremely useful tools if your language supports such things in a
reasonable manner. Poor implementations of these ideas shouldn't be
taken as proof that they aren't.

Bill Shortall

2004-12-13, 3:57 am

"Jentje Goslinga" <goslinga@telus.net> wrote in message
news:41BCA583.8050605@telus.net...
>
>
> Bill Shortall wrote:
Matrix[color=darkred]
and[color=darkred]
>
> This is not a very common operation which I have only
> encountered when working with block matrices.
> C++ is a powerful language but the simple algebra of
> Matrices, vectors and Scalars is sufficiently complex
> that using C++ shorthand notation is not always possible.
> Have you looked at OON, the Object Oriented Numerics
> Group?
>
> Akll of the basic Linear algebra
>
> You must be kidding, by the time I am done with Schur, SVD
> and Block Reflector, I will have 25-30k lines of C and a few
> thousand lines of Intel Assembler. And that is without Linear
> Equations or Complex Data type.
> And it is including testing modules for everything.
>
>
> Isn't that a bit large? My executable which contains all of
> the said factorizations and comprehensive testing software
> weighs in at 76 kilobytes. You generate 40 times the object
> code with one fifth the code, a factor 200!
>
play[color=darkred]
it[color=darkred]
>
> Have you tested those functions by constructing pathological
> matrices from their eigenvalues (for symmetric matrices)
> or from their singular values (in the case of the SVD) and
> likewise for Linear Equations and evaluated the accuracy?
>
>
> Jentje Goslinga

Subject: Re: Who uses clapack?

hi Jentje,

Bill wrote

class and then some more compilcated classes like a vector of vectors and

Jentje wrote
This is not a very common operation which I have only[color=darkred]
> encountered when working with block matrices.
> C++ is a powerful language but the simple algebra of
> Matrices, vectors and Scalars is sufficiently complex
> that using C++ shorthand notation is not always possible.

I believe that operator overloading is one of the primary advantages of
C++ The ability to write simple equations like A = B + C*D makes the
underlying mathematics stand out and the code easy to debug and remember.
I haven't had many problems wwith the shorthand notation. It is however
challenging to do in complex variables.

Bill wrote
ie solution of equations, decompositions, SVD eigenvectors etc. etc.

Jentje wrote
You must be kidding, by the time I am done with Schur, SVD[color=darkred]
> and Block Reflector, I will have 25-30k lines of C and a few

I am not kidding !!!
Remember I am writing templated code so that one function description
can in most cases handle float, double, complex<float>, and complex<double>
its all the same function....this reduces code length by a factor of 4.
Also remenber that I am using the standard template library which handles
the pointers and data types.. As am example

template <class T> inline
SMatrix<T>& SMatrix<T>::operator +=(const SMatrix<T> &B){
SYS_ASSERT(( (rows() == B.rows()) & (cols() == B.cols()) ),
"error += matrices not equal size \n");
std::transform(begin(),end(),B.begin(),begin(),std::plus<T>());
return *this;
}

here in 8 lines we have the code to do an update add += for all data types.
It also
checks for proper matrix size. The c code for this would be twice as long.
I also reuse the code in SVD's and eigenvectors to save space.
Bill
> and compiles into a static library of ~~ 3 megabytes

Jentje
> Isn't that a bit large? My executable which contains all of
> the said factorizations and comprehensive testing software

Remember that my lib file contains instances of code for all data types
For a particular problem only a small fraction of this is used and appears
in the object file of a particular app.

Jentje

Have you tested those functions by constructing pathological
matrices from their eigenvalues (for symmetric matrices)
or from their singular values (in the case of the SVD) and
likewise for Linear Equations and evaluated the accuracy?

Bill

I have tested all functions up to size 500 x 500 using
matrices filled with random numbers. I havent done anything
with pathological matrices yet.

I hope that this encourages others to consider using C++ for
their linear algebra problems

Regards...Bill Shortall

>

James Van Buskirk

2004-12-13, 3:57 am

"Jentje Goslinga" <goslinga@telus.net> wrote in message
news:41BD020D.8060503@telus.net...

> Yes, I had understood that, but don't you have to write
> separate implementations of those template functions for
> complex numbers for example. OK, I understand now, the
> complex float and complex double are the same code and the
> real float and double are another code, understood.

Now I don't think I get it. The Fortran version needs two
separate codes to handle the 6 possibilities: a real code
for single, double, and quad precision real and a complex
code for single, double, and quad precision complex. But
the previous poster seems to be saying that in the C++
version, a single code suffices for the four cases of
(single vs. double precision) and (real vs. complex).
There are at least some level 1 blas subprograms where
there are two complex versions and only one real version,
e.g. cdotc/zdotc/ydotc & cdotu/zdotu/ydotu for complex
data but only sdot/ddot/qdot for real data. It seems to
me that you aren't quite going to get a 4X code reduction
factor for several of the blas subprograms.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

kia

2004-12-13, 3:57 am

Victor Eijkhout wrote:
>
> The authors of lapack/scalapack are starting to work on a new release of
> these packages. One of the things we want to address is the sillyness
> that C users, instead of linking to the binaries of the Fortran version,
> use an automatically translated C version.

Your concerns seem oddly out of step, now that clc's own gatekeeper
declared them *off topic*. I agree, they should concentrate on doing
what C is good for -- numerical computing is not one of them.
Consequently, don't dilute your efforts to accommodate technically
illiterate crowd that could neither appreciate nor effectively use what
you have to offer.

Jan Vorbrüggen

2004-12-13, 8:57 am

> Interoperability between languages is a gray area where one
> needs to consider the following:
> 1 Ability to suppress C++ name decoration, knowledge of other
> conventions such as underscores appended or prepended to name
> 2 Knowledge of linker conventions such as upper case conversion
> and maximum name length
> 3 Ability for the programmer to specify parameter passing by
> address and by value which C and C++ have
> 4 Ability for the user to use parameter passing from left to
> right or vice versa
> 5 Conventions on stack based parameter passing notably padding
> and stack alignment
> 6 Decide on who resets the stack pointer, caller or callee.

You lack historical perspective. Coming from a DEC world, I was surprised
to learn that there existed systems where mixing source from different
languages in one application was not the norm, not say the rarest of
exceptions fraught with difficulties and pain - and all of that for no
perceptible reason at all but perhaps a bad case of NIH syndrome.

Although I must admit they didn't deem C++ important enough to re-write
the linker from scratch to avoid name mangling. However, the support tools
were quite good in that respect.

Jan

Jan Vorbrüggen

2004-12-13, 8:57 am

> Ok, it's clear where this use of temporaries comes from. However, it's
> not clear that dynamic allocation is the best solution. If your routine
> is called once, with a large dataset, then it could do its own
> alloation. However, if the routine is called a large number of times
> with the same small problem size, then you want to do the allocation
> outside it.

Just so. What I did in 187.facerec was to provide the storage by the caller
in form of an allocatable array*, and the called routined checked whether it
was 1) allocated and 2) of the proper size; if either of those conditions
was not met, it would (re-)allocate at that proper size.

Also, I believe some interfaces of that era - for instance, FFTPACK - used
the work arrays to store data structures describing the operation, or some
tables of constants, into them. These have more natural solutions in F95
that are much easier to use.

* I used POINTERs because most compilers didn't implement the ALLOCTABLE TR
at the time - I'm not even sure it was an official TR yet...Nowadays, used
ALLOCATABLE dummy arguments is the method of choice in almost all cases, IMO.

Jan

Jan Vorbrüggen

2004-12-13, 8:57 am

> In 2000 I had the misfortune to have to work Sun machines
> again, this time in C.
> The simulation model I was using gave quite different
> results on different machines with different versions
> of the libraries. I simply don't trust Sun floating point.

Then you should clearly get a new job, as you don't understand what
"numerical analysis" is about.

For the record, there's nothing wrong about FP on Sun machines, in general.
In some parts, they're even too close to the IEEE standard for predictable
performance.

Jan

Tino

2004-12-13, 4:08 pm

Bill Shortall wrote:
> "Victor Eijkhout" <see.sig@for.addy> wrote in message
> news:1gokv1s.16zyrkl1ytczkbN%see.sig@for.addy...
being[color=darkred]
LAPACK[color=darkred]
C++[color=darkred]
the[color=darkred]
will[color=darkred]
>
> Hi Victor,
>
> For the last few years I have been working on a
> set of C++ classes to do linear algebra. I was trying to
> acheive most of the functionality of Lapack/Blas but
> in a very user friendly fashion. The classes center
> arround a general purpose vector class along with a dense rectangular
Matrix
> class and then some more compilcated classes like a vector of vectors
and
> vector of matrices.
> All of them use operator overloading so you can write
> A = B + C * D where A,B,C, D are matrices and all of the
operators
> work for both real and complex matrices. Akll of the basic Linear
algebra
> functions are supported
> ie solution of equations, decompositions, SVD eigenvectors etc. etc.
> The entire library is less than 5000 lines of code
> and compiles into a static library of ~~ 3 megabytes
> You don't have to compile it but you will need either a
> Miccrosoft VC6 or Linux GCC compiler to use it.
> I have demo version available and if anyone in the group wants to
play
> with it and has a compatible compiler send me a line. Once you've
used it
> you will never go back to LAPACK
> pecos@cminet.net
> regards...Bill

I have also developed a C++ library which acts as nothing more than a
wrapper for the BLAS/LAPACK or Intel MKL libraries. In my case, I
claim only to support a subset of the BLAS/LAPACK functionality (though
it is 100% of the functionality that I use). I also support all of the
different matrix formats (dense, banded, symmetric), though only in
double precision (no complex matrices). The reason I was interested in
a more "official" effort is, that I believe that LAPACK's popularity
and utility is in its generality and I would welcome a sort of
standardization of the C++ interface as there seem to be dozens if not
hundereds of them floating around.

Ryan

Bill Shortall

2004-12-13, 4:08 pm

"Tino" <tino52@yahoo.com> wrote in message
news:1102947808.721378.166920@f14g2000cwb.googlegroups.com...
>
> Bill Shortall wrote:
> being
> LAPACK
> C++
> the
> will
> Matrix
> and
> operators
> algebra
> play
> used it
>
> I have also developed a C++ library which acts as nothing more than a
> wrapper for the BLAS/LAPACK or Intel MKL libraries. In my case, I
> claim only to support a subset of the BLAS/LAPACK functionality (though
> it is 100% of the functionality that I use). I also support all of the
> different matrix formats (dense, banded, symmetric), though only in
> double precision (no complex matrices). The reason I was interested in
> a more "official" effort is, that I believe that LAPACK's popularity
> and utility is in its generality and I would welcome a sort of
> standardization of the C++ interface as there seem to be dozens if not
> hundereds of them floating around.
>
> Ryan

Yes...the writers of LAPACK should write a C++
version. After all, they are the experts. It might even
be fun for them !
Bill ............. pecos@cminet.net
>

Victor Eijkhout

2004-12-13, 4:08 pm

Jan Vorbrüggen <jvorbrueggen-not@mediasec.de> wrote:

>
> Then you should clearly get a new job, as you don't understand what
> "numerical analysis" is about.

Harsh. Not all numerical algorithms automatically damp out
perturbations, and sometimes you just have to live with that.

V.
--
email: lastname at cs utk edu
homepage: www cs utk edu tilde lastname

Jan Vorbrüggen

2004-12-13, 4:08 pm

> Harsh. Not all numerical algorithms automatically damp out
> perturbations, and sometimes you just have to live with that.

Sure. But then put the blame on the hardware, instead of the
algorithm? That doesn't make any sense at all to me.

Jan

glen herrmannsfeldt

2004-12-13, 4:08 pm

(comp.lang.c removed)

kia wrote:

(snip)

> Your concerns seem oddly out of step, now that clc's own gatekeeper
> declared them *off topic*. I agree, they should concentrate on doing
> what C is good for -- numerical computing is not one of them.
> Consequently, don't dilute your efforts to accommodate technically
> illiterate crowd that could neither appreciate nor effectively use what
> you have to offer.

Topicality in comp.lang.c is very different than in some similar
groups, and almost has to be. Even so, there are too many posts
to read in a reasonable time. Note that the discussion of programs
written in Fortran is considered on topic in comp.lang.fortran, but the
discussion of C programs is not for comp.lang.c. That is, just being
written in C isn't enough. If the question is about the C language,
and even though it may require a sample program to understand the
question, then it is likely on topic.

Similarly, the comp.arch.powerpc newsgroup discusses the architecture
of the PowerPC processor, and not general questions related to macintosh
software as run on Apple machines using that processor. That doesn't
seem to stop people from posting them.

-- glen

Jentje Goslinga

2004-12-13, 4:08 pm

Jan Vorbrüggen wrote:
>
>
> Then you should clearly get a new job, as you don't understand what
> "numerical analysis" is about.
>
> For the record, there's nothing wrong about FP on Sun machines, in
general.
> In some parts, they're even too close to the IEEE standard for
predictable
> performance.

Your bitter and sarcastic attack serves no useful purpose.

This thread started with a request for suggestions for Lapack
and questions about the use of CLapack.
I have made some very useful suggestions for improvement of
Lapack which have been acknowledged by the OP.

I am developing C code for the most important decompositions
because I want them to be available in the C language in code
which can be easily maintained and also because I them in my
other work.
I have designed a recursive D&C solver for the Symmetric
Eigenvalue problem and am working on one for the SVD.
I have read several papers by Cuppen, Gu and Eisenstadt,
Tisseur and Dongarra, Gu, Demmel and Dhillon.
If I did not think I could write somewhat faster and somewhat
more accurate code than Lapack I would not be doing this.
Does this mean that I a "competitor" of Lapack?
No, I do this because I like to do it and will decide some
day what to do with it, and will probably end up giving it
away.

I really don't care whether some Jan Vorbruggen thinks that
I am neither a programmer nor a numerical analyst.

> Jan

Jentje Goslinga

Jentje Goslinga

2004-12-13, 4:08 pm

Victor Eijkhout wrote:
> Jan Vorbrüggen <jvorbrueggen-not@mediasec.de> wrote:
>
>
>
>
> Harsh. Not all numerical algorithms automatically damp out
> perturbations, and sometimes you just have to live with that.

This actually was a transient simulation model in C and it
does not take much for these simulations to diverge on
different platforms. On the other hand, I have tested complex
transient models between Intel and AMD processors and the
results were exactly the same.

Maybe the Sun platform is better now than at the time I used
it which is a long time ago.
Maybe the Lapack test suites should be expanded upon somewhat
and used to test the accuracy and speed of numerical linear
algebra between platforms. Using a set of difficult matrices
and a portable random generator.

> V.

Jentje Goslinga

2004-12-13, 9:01 pm

Victor Eijkhout wrote:
>
> However, it's
> not clear that dynamic allocation is the best solution. If your routine
> is called once, with a large dataset, then it could do its own
> alloation. However, if the routine is called a large number of times
> with the same small problem size, then you want to do the allocation
> outside it.

It isn't insofar as performance is concerned - a significant degradation
occurs just by introducing dynamic allocation. If that's acceptable then
the second part is a no brainer - always allocate inside and only if
size has changed - no need to treat large or small, frequent or glacial
calls separately.

Robert Corbett

2004-12-14, 3:57 am

In article <41BD36A3.4ED52EC6@sparrow.com>, kia <kia@sparrow.com> wrote:
>
>Your concerns seem oddly out of step, now that clc's own gatekeeper
>declared them *off topic*. I agree, they should concentrate on doing
>what C is good for -- numerical computing is not one of them.

While I might agree with you regarding C89, C90, or C95,
I disagree regarding C99. The facilities for numerical
computing in C99 are pretty much on a par with those in
Fortran, and in some cases are better. In particular,
C99 includes imaginary floating-point data types, which
are quite handy for some applications. It is still
easier to deal with arrays in Fortran than C, but C now
matches Fortran in areas where Fortran once had the edge,
such as aliasing and complex arithmetic.

Sincerely,
Bob Corbett

Sincerely,
Bob Corbett

Charles Russell

2004-12-14, 9:01 pm

"Robert Corbett" wrote

> It is still easier to deal with arrays in Fortran than C

and arrays of reals are what scientific computing is all about.

Jentje Goslinga

2004-12-14, 9:01 pm

Robert Corbett wrote:
> In article <41BD36A3.4ED52EC6@sparrow.com>, kia <kia@sparrow.com> wrote:
>
>
>
> While I might agree with you regarding C89, C90, or C95,
> I disagree regarding C99. The facilities for numerical
> computing in C99 are pretty much on a par with those in
> Fortran, and in some cases are better. In particular,
> C99 includes imaginary floating-point data types, which
> are quite handy for some applications. It is still
> easier to deal with arrays in Fortran than C, but C now

<snip>

> Sincerely,
> Bob Corbett

I agree with what you are saying and have a few more notes
on t