Newsgroup sci.math.num-analysis 28902

Directory

dlundin@aol.com wrote:
: Gentlemen,
: I am attempting to use the SVD code from www.netlib.org for CLAPACK, the C
: language version of LAPACK.
: I have discovered that several calls are made to routines which I do not
: have and did not find any directions that would suggest I need to download
: any other libraries.
: The routines (so far) are:
: dcopy, dgemm_, dgemv_, dlamc1_, dlamc3_, dlamc4_, dlamc5_, dlarfg_,
: dlarf_, dlasq2_, dlasrt_, dlam2r_, dorml2_, drot_, dscal_, dswap_, dtrmm_,
: dtrmv_, d_sign, pow_dd, pow_di, s_cat, s_cmp, s_copy.
: Are there any gurus out there who can assist?
looks like BLAS-routines....
uwe.
====================================================================
=== "the practical scientist is trying to solve tomorrow's       ===
===  problem with yesterday's computer;                          ===
===  the computer scientist often has it the other way around."  ===
====================================================================
===       Uwe Schmitt, Bahnhofstr. 79, 66111 Saarbruecken        ===
===       phone 0681 / 37 36 35                                  ===
===       email schmitt@num.uni-sb.de                            ===
====================================================================

In article <3295C259.26DB@oitc.com>, Tom Shaw   wrote:
>Hans D. Mittelmann wrote:
> Also don't forget that fortran and c handle multidimensional arrays
> differently.
And C even handles multi-dimensional arrays different from C....
DId this sound cryptic?  That's because C has two different ways
to simulate multi-dimensional arrays .... the language itself
doesn't have them built-in really.
Even single-dimensional arrays are "second-class citizens" in C
-- try this example:
    void sub( int arr[10] )
    {
        int i;
        for( i = 0; i < sizeof(arr)/sizeof(arr[0]); i++ )
            printf( " %d", arr[i] );
    }
    main()
    {
        int arr[10] = { 1,2,3,4,5,6,7,8,9,10 };
        sub( arr );
    }
This program does not output  " 1 2 3 4 5 6 7 8 9 10" as one naively
might expect.  Instead it outputs " 1" or " 1 2", depending on the
environment used.....
-- 
----------------------------------------------------------------
Paul Schlyter,  Swedish Amateur Astronomer's Society (SAAF)
Grev Turegatan 40,  S-114 38 Stockholm,  SWEDEN
e-mail:  pausch@saaf.se     psr@home.ausys.se    paul@inorbit.com

In article <57ge2s$58r@news.ox.ac.uk>, trin0008@ermine.ox.ac.uk (Rick Heylen) writes:
|> I've implemented the routines tred2 and tqli in Java. These routines 
|> tridiagonalise a symmetric matrix and diagonalise a tridiagonal one 
|> respectively. The book says that "Complex versions of the previous 
|> routines .. tred2, and tqli are quite analogous to their real counterparts."
|>           .    .    .
|> Now I'm trying to avoid having to go to the bother of understanding these 
|> crappily written routines because that's the general idea of providing them.
|>           .    .    .
|> Please do not direct me to CLAPACK as porting a machine's translation 
|> into C of hand optimised FORTRAN 70 into Java would drive me round the bend.
|> It's bad enough having the Numerical Recipes arrays starting at 1 and the 
|> code liberally sprinkled with breaks continues and massively block 
|> structures ifs and the impenetrable variable naming and .... etc 
Well, if you insist on starting with junk and try to build serious
software, you will either fail or go through hell.  You brought your
problems on yourself.  Hand your copy of Numerical Recipes over to
the Public Hangman (who was the official responsible for burning
obscene publications) and start from a better basis.
Try looking at TOMS, or even go back to the original Wilkinson and
Reinsch work.  And why did you ignore LAPACK itself, which is in
fairly clean Fortran 77 and well (if somewhat BLASarrly) structured.
One could ask why you are trying to perform complex matrix algebra
in Java, but I shall desist :-)
Nick Maclaren,
University of Cambridge Computer Laboratory,
New Museums Site, Pembroke Street, Cambridge CB2 3QG, England.
Email:  nmm1@cam.ac.uk
Tel.:  +44 1223 334761    Fax:  +44 1223 334679

Au_Yeung_Man_Ching wrote:
> 
> Hi,
> 
> I would like to ask some questions:
> 
> (1) Given the following problem,
> 
>        min J(q,S)
>        q,S
> 
>     in doing the above problem, we do the
>     above problem like this:
> 
>      min min J(q,S)
>       q   S
> 
>     i.e. doing the problem sequentially as two
>     parallel independent minimization problem.
>     With a fixed q first, minimize J(q,S) w.r.t.
>     S to obtain S*. Then, minimize J(q,S*) w.r.t.
>     q. Finally, get the optimal point(local/global??)
>     (q*,S*).
> 
>     CAN WE DO THE PROBLEM AS DESCRIBED ABOVE??
>     IF CANNOT, how can we do?
> 
> Thanks in advance!!
> Regards,
> --
> +------------------------+--------------------+
> | Au Yeung Man Ching     | Rm 304             |
> | Electronic Engineering | Ho Sin Hang Bldg.  |
> | CUHK                   | CUHK               |
> +------------------------+--------------------+
> | Email : mcau@ee.cuhk.edu.hk                 |
> | Homepage : http://www.ee.cuhk.edu.hk/~mcau/ |
> +---------------------------------------------+
> 
> --
> +------------------------+--------------------+
> | Au Yeung Man Ching     | Rm 304             |
> | Electronic Engineering | Ho Sin Hang Bldg.  |
> | CUHK                   | CUHK               |
> +------------------------+--------------------+
> | Email : mcau@ee.cuhk.edu.hk                 |
> | Homepage : http://www.ee.cuhk.edu.hk/~mcau/ |
> +---------------------------------------------+
Hi,
to solve a problem by alternatingly minimizing with respect to one of
the parameters is called "cyclic coordinate descent". It converges
linearly under standard hypotheses. See, for example, D.G. Luenberger,
Linear and Nonlinear Programming, p 227-230.
Hope that helps
-- 
Hans D. Mittelmann			http://plato.la.asu.edu/
Arizona State University		Phone: (602) 965-6595
Department of Mathematics		Fax:   (602) 965-0461
Tempe, AZ 85287-1804			email: mittelmann@asu.edu

Ron Shepard wrote:
> 
> In article <3299C149.3A9659B2@asu.edu>, Hans D Mittelmann
>  wrote:
> 
> [...]
> >here is my theorem. Unless it can be found in the literature, I'd like
> >to be quoted as the source.
> >
> >    D det(a) = sum(i=1,n) product(j.ne.i) lambda_j(A)
> >
> >This requires one call of, say, the QR algorithm and is thus a O(n^3)
> >method.
> 
> I think there is a factor missing.  I get:
> 
>     D det(A) = sum(i=1,n) (D(lambda_i) Product(j.ne.i) lambda_j )
> 
> and therefore you need the derivatives of the eigenvalues, which can be
> computed from the perturbation-theory expression (I'm assuming all the
> derivatives are evaluated at the same reference point, 0; oh yeah, and the
> matrix A is normal).  I'm not sure this is the best way of solving this
> problem, but the above expression points out that the value is zero if
> there is more than one zero eigenvalue (the product terms will always be
> zero in this case, having at least 1 and sometimes more zero factors), and
> the expression is zero if the zero eigenvalue is quadratic (i.e. O(2) or
> higher, the D(lambda) term is zero for the surviving product term in this
> case; BTW, I think this must be true if A is semidefinite or if its
> inertia triplet is locally constant).  So in case both of these conditions
> are satisfied, there is only one term that contributes to the summation,
> 
>    D det(A) = D(lambda_i) Product(j.ne.i) lambda_j  ;where lambda_i=0
> 
> which is why I suspect that this may not be the best computational
> approach (too much effort, too simple a result ;-).
> 
> $.02 -Ron Shepard
Hi,
I think, Ron, you believe I was just differentiating the expression
that  the det is the product of the eigenvalues. I absolutely was not!
Note that the original question was not directly for the derivative of
det(A) but for a more economical way to compute trace (B det(A)). That
is the expression I rewrote (with one little trick). So, my result holds
for any rank of A and it is NOT true that D det(A) is zero as soon as A
has at least one zero eigenvalue.
Hans
-- 
Hans D. Mittelmann			http://plato.la.asu.edu/
Arizona State University		Phone: (602) 965-6595
Department of Mathematics		Fax:   (602) 965-0461
Tempe, AZ 85287-1804			email: mittelmann@asu.edu

> Cedric Dourthe wrote:
>   
> > I'm looking for iterativ methods package in order to solve linear problems
> > I found C and fortran routines but only for sparse matrix
> > 1) where can I find the same packages for dense and complex matrix?
> > 2) Comparisons iterative methods/direct methods for this type of matrix?
  The Netlib repository contains the subroutine library QMRPACK,
  which contains various versions of the QMR iterative method for
  solving linear systems. There are subroutines also for complex 
  matrices.
  In these subroutines the matrix-vector multiplication is external 
  to the solvers. This means that the algorithm routines will return to the 
  caller each time a matrix-vector multiplication is required, signal to 
  the caller by means of a flag whether they need a multiplication by A or 
  by A^T (where appropriate), and then expect the caller to call them back 
  after the multiplication is done.
  So, you may have either a sparse or a dense matrix.
  The Netlib is available from  http://www.netlib.org/
*... Mr. Erkki M. S. Heikkola      emsh@math.jyu.fi .............*
*... Department of Mathematics     http://www.math.jyu.fi/~emsh .*
*... University of Jyvaskyla       tel. +358 (0)14 645651 .......*
*... P.O. Box 35, 40351 Jyvaskyla       +358 (0)14 602760 .......*
*... FINLAND                       fax  +358 (0)14 602731 .......*

In article <329C4896.F0DA37D@asu.edu>,
Hans D Mittelmann   wrote:
>Ron Shepard wrote:
>> 
>> I think there is a factor missing.  I get:
>> 
>>     D det(A) = sum(i=1,n) (D(lambda_i) Product(j.ne.i) lambda_j )
>> and therefore you need the derivatives of the eigenvalues, which can be
>> computed from the perturbation-theory expression (I'm assuming all the
>> derivatives are evaluated at the same reference point, 0; oh yeah, and the
>> matrix A is normal).  I'm not sure this is the best way of solving this
>> problem, but the above expression points out that the value is zero if
>> there is more than one zero eigenvalue (the product terms will always be
>> zero in this case, having at least 1 and sometimes more zero factors), and
>> the expression is zero if the zero eigenvalue is quadratic (i.e. O(2) or
>> higher, the D(lambda) term is zero for the surviving product term in this
>> case; BTW, I think this must be true if A is semidefinite or if its
>> inertia triplet is locally constant).  So in case both of these conditions
>> are satisfied, there is only one term that contributes to the summation,
>> 
>>    D det(A) = D(lambda_i) Product(j.ne.i) lambda_j  ;where lambda_i=0
>> 
>> which is why I suspect that this may not be the best computational
>> approach (too much effort, too simple a result ;-).
>> 
>> $.02 -Ron Shepard
Note that given an eigenvalue decomposition of A from the QR-algorithm
the perturbation analysis seems to be really simple.
Assume that A is diagonalizable and X^(-1) A X has the eigenvalues on
the diagonal, and X^(-1) DA X=d_{i,j}.
The result is then: 
  D det(A)=sum(i=1,n) d_{i,i}* product(j.ne.i) lambda_j(A)
Computable, but expensive (you also have to drop the imaginary part
which in in-exact arithmetic will be slightly non-zero).
(And as you say if lamba_i=0 it gets simpler).
>Hi,
>I think, Ron, you believe I was just differentiating the expression
>that  the det is the product of the eigenvalues. I absolutely was not!
Your result indicate that you did.
>Note that the original question was not directly for the derivative of
>det(A) but for a more economical way to compute trace (B det(A)).
The original question was about D det A=trace(B.DA)
I assume that D stands for differential and it is thus a question about
evaluating trace(B.DA), which obviously must should contain DA.
>That
>is the expression I rewrote (with one little trick). So, my result holds
>for any rank of A and it is NOT true that D det(A) is zero as soon as A
>has at least one zero eigenvalue.
Ron said it was zero if A had more than one zero eigenvalue.
--
// Homepage  http://www.dna.lth.se/home/Hans_Olsson/
// Email To..Hans.Olsson@dna.lth.se [Please no junk e-mail]

Harold Howe wrote:
> 
> I need a peak finding algorithm, and I have no idea where to
> start looking for one.  What are some of the techniques used
> to analyze sampled data.
> 
> If you have any advice I would appreciate it.
> 
I think the beasts you seek are the Fibonacci and Golden
Section Search methods.  See D G Luenberger, "Linear
and Non-Linear Programming".
I also use a method where I examine the function at three 
equally spaced points to give y1 = y(x1), and so forth
for y2 and y3.  Then use some logic of the sort:
   if (y1 < y2 and y2 < y3) {
       y1 = y2;
       y2 = y3;
       y3 gets bigger;
   }
and so forth.  
You have to be able to slide all 3 points to the left or right,
as above, and also contract the points if the middle one
is the largest of the three.
The referenced methods are similar, but faster and more elegant
(read "tedious" by the person who does the coding).
Selecting the initial search region can be troublesome...
-- 
***************************************************************
  Dr. Albert P. Gerheim   |    Sonalysts, Inc. PO. Box 280 
  gerheim@sonalysts.com   |         215 Parkway North
  1 (800) 526-8091 X 218  |        Waterford, CT  06385
***************************************************************

Hans Olsson wrote:
> 
> In article <329C4896.F0DA37D@asu.edu>,
> Hans D Mittelmann   wrote:
> >Ron Shepard wrote:
> >>
> >> I think there is a factor missing.  I get:
> >>
> >>     D det(A) = sum(i=1,n) (D(lambda_i) Product(j.ne.i) lambda_j )
> >> and therefore you need the derivatives of the eigenvalues, which can be
> >> computed from the perturbation-theory expression (I'm assuming all the
> >> derivatives are evaluated at the same reference point, 0; oh yeah, and the
> >> matrix A is normal).  I'm not sure this is the best way of solving this
> >> problem, but the above expression points out that the value is zero if
> >> there is more than one zero eigenvalue (the product terms will always be
> >> zero in this case, having at least 1 and sometimes more zero factors), and
> >> the expression is zero if the zero eigenvalue is quadratic (i.e. O(2) or
> >> higher, the D(lambda) term is zero for the surviving product term in this
> >> case; BTW, I think this must be true if A is semidefinite or if its
> >> inertia triplet is locally constant).  So in case both of these conditions
> >> are satisfied, there is only one term that contributes to the summation,
> >>
> >>    D det(A) = D(lambda_i) Product(j.ne.i) lambda_j  ;where lambda_i=0
> >>
> >> which is why I suspect that this may not be the best computational
> >> approach (too much effort, too simple a result ;-).
> >>
> >> $.02 -Ron Shepard
> 
> Note that given an eigenvalue decomposition of A from the QR-algorithm
> the perturbation analysis seems to be really simple.
> Assume that A is diagonalizable and X^(-1) A X has the eigenvalues on
> the diagonal, and X^(-1) DA X=d_{i,j}.
> The result is then:
> 
>   D det(A)=sum(i=1,n) d_{i,i}* product(j.ne.i) lambda_j(A)
> 
> Computable, but expensive (you also have to drop the imaginary part
> which in in-exact arithmetic will be slightly non-zero).
> (And as you say if lamba_i=0 it gets simpler).
> 
> >Hi,
> >I think, Ron, you believe I was just differentiating the expression
> >that  the det is the product of the eigenvalues. I absolutely was not!
> 
> Your result indicate that you did.
> 
> >Note that the original question was not directly for the derivative of
> >det(A) but for a more economical way to compute trace (B det(A)).
> 
> The original question was about D det A=trace(B.DA)
> 
> I assume that D stands for differential and it is thus a question about
> evaluating trace(B.DA), which obviously must should contain DA.
> 
> >That
> >is the expression I rewrote (with one little trick). So, my result holds
> >for any rank of A and it is NOT true that D det(A) is zero as soon as A
> >has at least one zero eigenvalue.
> 
> Ron said it was zero if A had more than one zero eigenvalue.
> 
> --
> // Homepage  http://www.dna.lth.se/home/Hans_Olsson/
> // Email To..Hans.Olsson@dna.lth.se [Please no junk e-mail]
Hi,
I just noticed my error before I saw your message. The notation "D"
simply confused me and I was expressing trace (B detA) instead of trace
(B DA). The questioner may have a solution at least for a deficiency of
exactly 1. (I'll be gone for a week now).
Hans 
-- 
Hans D. Mittelmann			http://plato.la.asu.edu/
Arizona State University		Phone: (602) 965-6595
Department of Mathematics		Fax:   (602) 965-0461
Tempe, AZ 85287-1804			email: mittelmann@asu.edu

SOME NEW RESULTS IN THE FIELD OF DISCRETE MATH/DESIGNS/CODES
Date  : 96/11/27 We
Author: Uenal Mutlu (bm373592@muenchen.org)
The following bounds of enumerated 'm=t COVERING DESIGNS' are to
the best of my knowledge *NEW* upper bound results compared to
known published or online available results (cf. [1],[2],[4]).
Definition:
 A Covering Design C(v,k,t,m,l,b) is a pair (V,B), where V is a set
 of v elements (called points) and B is a collection of b k-subsets
 of V (called blocks), such that every m-subset of V intersects at
 least l members of B in at least t points (v >= k >= t and m >= t).
New upper bounds for C(v,k,t,m,l=1,b). Listed are m=t designs only:
 k t m	v     b (bOld)	Remarks/Author/Method
------------------------------------------------------------------
 6 4 4 15   118 (120)	Rade Belic, [3]
 6 5 5 14   377 (378)	Rade Belic, [3]
 6 5 5 15   609 (610)	Combine-Construction [1]
 7 5 5 17   408 (463)	Rade Belic, [3]
 7 5 5 18   618 (663)	Combine-Construction [1]
 7 5 5 19   772 (824)	Turan Theory/Sidorenko Construction [1]
 7 5 5 20  1115 (1153)	Combine-Construction [1]
 7 5 5 21  1424 (1476)	[3]
 7 5 5 22  1938 (1990)	[3]
 7 5 5 23  2405 (2458)	[3]
 7 5 5 24  3096 (3208)	[3]
 7 5 5 25  3772 (3830)	[3]
(bOld refers to the upper bounds of [2] on 96/11/26)
The new upper bounds were obtained by some published (see [1]) and
some not yet published algorithms and programs of Rade Belic [3]
and Uenal Mutlu (including implementations of some recursive
constructions from [1]).
These results will also improve some other designs with higher
k,t,m if at least some of the well known recursive construction
methods are applied.
Not shown are the latest improvements and new upper bounds for
designs with m > t. A list containing all upper bounds of enumerated
designs in the range k <= 7, v <= 54, t <= 7, m <= 7 and b <= 9999
can be requested from the author [5]).
Pointers to further improvements and new results are welcome.
Uenal Mutlu
References:
[1] D.M.Gordon, O.Patashnik, G.Kuperberg "New Constructions for
    Covering Designs", J.Combin.Designs 3.4 (1995) 269-284
    (for updates and errata see [2])
[2] Dan Gordon's "La Jolla Covering Repository" at
    http://sdcc12.ucsd.edu/~xm3dg/cover.html
    (as of 96/11/26)
[3] SOPT v1.0 07/96 (96/11/19) - Design Optimization Program
    written by Rade Belic
[4] D.R.Stinson "Coverings" in C.J.Colbourn, J.H.Dinitz (eds.)
    "The CRC Handbook of Combinatorics" CRC Press (1996) 260-265
[5] Uenal Mutlu "List of Covering Designs" (enumerated upper bounds)
    Includes also m > t designs. Current list is LST1127A.ZIP
-- Uenal Mutlu (bm373592@muenchen.org)   
   Math Research, Designs/Codes, Data Compression Algorithms, C/C++
   Loc  : Istanbul/Turkey + Munich/Germany

In article <329673b2.914629444@news.entertain.com>, davismic@entertain.com
(Michael Davis) wrote:
> I am looking for a shareware program for Windows 3.11 or DOS to
> calculate distances between two positions on the earth surface and the
> required azimuth to get there.  Does anyone know of one and where I
> can find it.  Please reply to the newsgroup or by Email to
> "davismic@entertain.com".  Thank you,  Michael Davis
ftp://ftp.netcom.com/pub/hb/hbaker/FAQ-lat-long.txt

RESERVOIR SIMULATION SOFTWARE DEVELOPMENT
Scientific Software-Intercomp, Inc. a worldwide leader in
software development and technical consulting for the petroleum
industry has an immediate opening for a Software Development
Engineer.
Candidate will be a recent graduate or soon-to-be graduate to
help with the development and maintenance of our reservoir
simulation software.  Individual must have skills in Fortran
programming, numerical analysis and  physics, chemistry or
engineering.  Coursework or experience in reservoir simulation
development, computational fluid dynamics and C or C++
programming would be a plus.   Individual will have a degree in
Computer Science, Math, Physics, Chemistry, or Engineering.
SSI offers a competitive salary and an outstanding employee
benefits package.  Please mail or fax your resume to:
                 Vice President, Human Resources
              Scientific Software- Intercomp, Inc.
                  1801 California St., Ste. 295
                        Denver, CO  80202
                      phone: (303) 292-1111
                       fax: (303) 295-2235
        SSI IS PROUD TO BE AN EQUAL OPPORTUNITY EMPLOYER

Hello!
I am looking for literature on adaptive grid generation. Both
static and dynamic. I am a beginer, so I would appreciate,
relatively easier materials first. Any references for books
and journal articles will be greatly appreciated.
Thank you. Regards,
		Nash'at
......................................................................
Nash'at Ahmad			|  O Sternenfall,
Desert Research Institute	|  von einer Brucke einmal eingesehn-: 
nashata@sage.dri.edu		|  Dich nicht vergessen. Stehn!
+(702) 324-1635			|             - Rainer Maria Rilke
......................................................................

An important point is the expected error in the original data. If
the error is very small, then you want to do interpolation. If the
error is substantial (e.g. 5% or more) then you may want to try
smoothing. Since you say the data is samples from some real world
process, I assume the second case is what you want.
Since the data is 2 dimensional, I would start by making a plot.
If the underlying dependence is "simple", then a least squares
polynomial fit would be the easiest thing to do. Numerical Recipes
(referenced in another message) has code to do that.
If the dependence is not simple, or very non-linear, then you may
first want to try transforming the independent or dependent variables
in order to make it simple. For intensity-like data, a log transform
on the dependent variable may be appropriate.
If the data resists all attempts at simplification, you may want
to resort to "smoothing spline" approaches, e.g. loess. You can find it
at:
	http://www.netlib.org/a/index.html
I hope this is helpful,
Sani.
-- 
----------------------------------------------------------------------
Sani R. Nassif         nassif@austin.ibm.com    EEEE EEEE  E   E
IBM Austin Research Lab                          EE  E   E EE EE
11400 Burnet Rd.  MS/9460                        EE  EEEE  E E E
Austin, TX 78758-3493  tie.      678-8161        EE  E   E E   E
tel. (512)838-8161     fax. (512)838-4036       EEEE EEEE  E   E
----------------------------------------------------------------------

Does anyone know where I might obtain the source code for an
implementation of the Fletcher-Xu method for solving nonlinear
least squares problems?  (as in Fletcher's "Practical Methods of
Optimization", 2nd Edition)  I would prefer something in Pascal,
but C/C++ or Fortran would certainly be okay too.
Thanks for any and all pointers on this subject.
Mark Vaughan
m.a.vaughan@larc.nasa.gov

In article <576rfd$lqg@electra.saaf.se>, pausch@electra.saaf.se says...
> 
>Even single-dimensional arrays are "second-class citizens" in C
>-- try this example:
> 
>    void sub( int arr[10] )
>    {
>        int i;
>        for( i = 0; i < sizeof(arr)/sizeof(arr[0]); i++ )
>            printf( " %d", arr[i] );
>    }
> 
>    main()
>    {
>        int arr[10] = { 1,2,3,4,5,6,7,8,9,10 };
>        sub( arr );
>    }
> 
>This program does not output  " 1 2 3 4 5 6 7 8 9 10" as one naively
>might expect.  Instead it outputs " 1" or " 1 2", depending on the
>environment used.....
> 
>-- 
>----------------------------------------------------------------
>Paul Schlyter,  Swedish Amateur Astronomer's Society (SAAF)
>Grev Turegatan 40,  S-114 38 Stockholm,  SWEDEN
>e-mail:  pausch@saaf.se     psr@home.ausys.se    paul@inorbit.com
The reason why you are getting the results is that 
sizeof (arr)/sizeof (arr[0]) = 1 not 10 as you might expect.
Try instead i<10, You will get the output "1 2 3 4 5 6 7 8 8 10".
Hyun M. Peck, Ph.D., Mechanical Engineer.  
The University of MIchigan

>>So, my result holds
>>for any rank of A and it is NOT true that D det(A) is zero as soon as A
>>has at least one zero eigenvalue.
>
>Ron said it was zero if A had more than one zero eigenvalue.
I was in a bit of a hurry when I posted originally so I should have been
clearer about this.  A _must_ have one zero eigenvalue in order to be
singular.  I pointed out, correctly I think, that the summation in the
general expression for D det(A) reduces to a single term when A has a zero
eigenvalue.  That single term also vanishes in some cases, such as when
there are additional zero eigenvalues, or when the original zero
eigenvalue is stable.
There are other ways to approach this problem too.  It seems to me that
the final expression seems too "simple" compared to the preceeding effort
and there may be other ways to compute the result with less effort.  For
example, the identity
    det( exp(B) ) = exp( Tr(B) )
or one of several like it, might be used somehow to arrive at a simpler
expression for the special case when A is singular. 
$.02 -Ron Shepard

Hyun Min Peck wrote:
> 
> In article <576rfd$lqg@electra.saaf.se>, pausch@electra.saaf.se says...
> >    void sub( int arr[10] )
> >    {
> >        int i;
> >        for( i = 0; i < sizeof(arr)/sizeof(arr[0]); i++ )
> >            printf( " %d", arr[i] );
> >    }
> >
> >    main()
> >    {
> >        int arr[10] = { 1,2,3,4,5,6,7,8,9,10 };
> >        sub( arr );
> >    }
> >
> >This program does not output  " 1 2 3 4 5 6 7 8 9 10" as one naively
> >might expect.  Instead it outputs " 1" or " 1 2", depending on the
> >environment used.....
> >
> The reason why you are getting the results is that
> sizeof (arr)/sizeof (arr[0]) = 1 not 10 as you might expect.
> Try instead i<10, You will get the output "1 2 3 4 5 6 7 8 8 10".
> 
Yes, of course that's what is happening. But the point is you
would not expect that to happen -- at least I didn't.
Add the following in the main routine and see what you get:
	printf("The size is %d\n",sizeof(arr)/sizeof(arr[0]));
I get 10. 
The same line added to the subroutine gives prints 1.
What's happening of course is that C is printing the size of 
the pointer passed to sub rather than the size of the array
that is declared as being passed.
For the record, I'm using the PC version of gcc (DJGPP).
---
Jim

Greetings!
I am trying to understand the origin and the rationale behind
some code that implements the following:
There is a comment there that says 
        Approximate Bias Correction
There are also comments that say "calculate throw".
Here is what the algebraic part of the code:
GIven:
        r = correlation
        error_threshold
        error_in_correlation
Calculated:
        N = (error_threshold * error_in_correlation)
        D = ( 1 - r** 2 ) = ( 1 -r) * ( 1 + r)
        E = N / D
        Z =(0.5*log((1+r)*(1-r)))+(RandomNUmber*E)+(r*E*E*0.5)
        R = ( Z -1 ) / (Z+1)
Could you help me with some pointers for this sort of correction
methods. Are there references out there that I could lookup.
yahoo and friends are too much to handle ... they give me everything
in the world. Any special email lists that I can subscribe to
better understand this. Thank you.
Please email to me directly.. I may not have access to this news group
for the next several days due to holidays.
Thanks and I greatly appreciate your time and generosity.
--
Raman
-------------------==== Posted via Deja News ====-----------------------
      http://www.dejanews.com/     Search, Read, Post to Usenet

I have a number of x,y data points that form an arc of constant radius. 
 Do you know a good technique for finding the center and radius of the 
arc?
I've tried the following, but I'm not sure of it's reliability:
take first data point, find line running through point but also tangent 
to the arc, take second data point and find line tangent again.  
Normalize both lines and find the intersection.  This will give me the 
center, from here the radius is easy.  
Do you have any suggestions?  
Thanks, Marc

On Thu, 14 Nov 1996, Channing Walton wrote:
> Hi,
> I developed a method of of integration which seems to be very quick and accurate and
> was wondering what problems may occur with it or if it has been done before (very 
> likely).
> It goes like this:
> You wish to integrate a function between A and B,
> fit a polynomial (5th order seems good) to the function and  you have your first
> approximation.
> then repeat by integrating between A and (B-A)/2 + integral between (B-A)/2 to B
> if the result is within a suitable tolerance to the first attempt then end
> if it is not then continue dividing and checking.
> 
> the point been that incoding this you can use a simple recursive procedure.
> 
> Now there are obvious problems with particular functions but if you ensure that the
> integral is chopped up into sufficient pieces then all is well.
> 
> I have found this to be very quick.
> 
> Channing
This isn't anything new, and there are many other quicker and more
accurate numerical integration methods.  
Some problems - polynomial interpolation of the function is a very bad way
of approximating the integrand.  Also, this is only good for nice
functions.
You may wish to look up other advanced methods of integration, such as
Adaptive Simpson's method or the method of Gaussian Quadrature.
Good luck.

Jim Hunter  writes:

> > >    void sub( int arr[10] )
> > >    {
> > >        int i;
> > >        for( i = 0; i < sizeof(arr)/sizeof(arr[0]); i++ )
> > >            printf( " %d", arr[i] );
> > >    }

> What's happening of course is that C is printing the size of 
> the pointer passed to sub rather than the size of the array
> that is declared as being passed.
Correct concept, but wrong wording: you can never define the type of a
parameter as an array. There is a catch-all rule which says that if
the type of a parameter as declared turns out to be an array type, its
type is changed to be a pointer to the element type instead. So, the
type of arr above is not `int[10]', but rather `int*'.
The same rule applies for parameters which seem to be declared as
functions: there actual type is a pointer to function instead.
Cheers
Tanmoy
-- 
tanmoy@qcd.lanl.gov(128.165.23.46) DECNET: BETA::"tanmoy@lanl.gov"(1.218=1242)
Tanmoy Bhattacharya O:T-8(MS B285)LANL,NM87545 H:#9,3000,Trinity Drive,NM87544
Others see ,
or. -- 
fax: 1 (505) 665 3003   voice: 1 (505) 665 4733    [ Home: 1 (505) 662 5596 ]

Hans D Mittelmann (mittelmann@asu.edu) wrote:
: Jenny Sanchez wrote:
: > 
: > i just downloaded dgels.c from netlib. I got the version that supposedly
: > includes all dependencies. I also got the f2clibs.
: > 
: > however the routines are still calling subroutines i didn't get.
: > These are
: > [...]
: > 
: > I've searched around the clapack offerings, and can't find these. anyone
: > know where these routines are?
: > 
: > thanks, Jenny
: Hi,
: these are BLAS routines. Get the BLAS lib, it's in CLAPACK also.
: 
: Hans
: -- 
Actually, you should probably first try to link to the optimized BLAS in your
system (link with -lblas), which is usually provided by the vendor.  It
makes a significant difference in terms of speed if you use the BLAS version 
optimized for your machine.  Of course, if you don't have it, you'll have
to download the routines from NETLIB and recompile them.
-Madhu
-------------------------------------------------------------------------
Madhu Nayakkankuppam                              Phone: 212 - 998 - 3350
Courant Institute of Mathematical Sciences        Fax: 212 - 995 - 4121
251 Mercer Street                                 Email: madhu@cs.nyu.edu
New York, NY 10012-1185     Web: http://www.cs.nyu.edu/phd_students/madhu
-------------------------------------------------------------------------

Hans D Mittelmann (mittelmann@asu.edu) wrote:
: Jenny Sanchez wrote:
: > 
: > i just downloaded dgels.c from netlib. I got the version that supposedly
: > includes all dependencies. I also got the f2clibs.
: > 
: > however the routines are still calling subroutines i didn't get.
: > These are
: > [...]
: > 
: > I've searched around the clapack offerings, and can't find these. anyone
: > know where these routines are?
: > 
: > thanks, Jenny
: Hi,
: these are BLAS routines. Get the BLAS lib, it's in CLAPACK also.
: 
: Hans
: -- 
Actually, you should probably first try to link to the optimized BLAS in your
system (link with -lblas), which is usually provided by the vendor.  It
makes a significant difference in terms of speed if you use the BLAS version 
optimized for your machine.  Of course, if you don't have it, you'll have
to download the routines from NETLIB and recompile them.
-Madhu
-------------------------------------------------------------------------
Madhu Nayakkankuppam                              Phone: 212 - 998 - 3350
Courant Institute of Mathematical Sciences        Fax: 212 - 995 - 4121
251 Mercer Street                                 Email: madhu@cs.nyu.edu
New York, NY 10012-1185     Web: http://www.cs.nyu.edu/phd_students/madhu
-------------------------------------------------------------------------

Has there been any good methods developed for numerical solution of
DAE's?  Is there anywhere on the net that has information about it?  I
read a book recently and it said there wasn't a good method for this
yet.  The book is a few years old, so I'm wondering what the current
situation is.
Thanks,
-joel
-- 
taotree Tutor and Stuff
Math and Physics Solver and thoughts on Creativity
http://www.geocities.com/CapeCanaveral/8103/

In article <5731q0$pkb@harbinger.cc.monash.edu.au>,
suter@fawlty8.eng.monash.edu.au (Mr D. Suter) wrote:
% In article  aroberts@usq.edu.au
(Tony Roberts) writes:
% >In article <3287C6B6.224@sumppu.fmi.fi>, Olaf Amm  wrote:
% >
% >> Hi,
% >> 
% >> I am looking for an algorithm that is able to estimate the curl and
% >> divergence of a vector field from measurements of that field, given on
% >> irregularly spaced points.
% >> It would be best if that algorithm could do it without gridding the data
% >
% >Could do a radial basis function fit of each component, then differentiate
% >the fits.  
% >
% >For scalar data u_i at n points x_i (in m-dimensions) fit
% >u(x)=\sum_{j=1}^n a_jf(|x-x_j|)  +b.x +c
% >where f(r)=r^2\log r or f(r)=r^3.
% >
% >There are n+m+1 unknowns: a_j, vector b, c:
% >n equations come from u(x_i)=u_i
% >m equations from \sum_{j=1}^na_jx_j=0
% >and 1 from   \sum_{j=1}^n a_j=0
% >This is a symmetric system of equations.  For economy chose about a dozen
% >data points in the neighbourhood of the place you want to estimate the
% >derivatives.
% 
% If you want to go this way, may I suggest that you look at splines
% that may "adapt" to the amount of div and curl (or the variations in these)
% in your data. A starting point in the Div-Curl splines of Amodei and
Benbourhim
% (J Approx Theory 1991). The get yourself a spline solving package
% like Chong Gu's rkpk. I'm not necessarily suggesting this will always
% be the best way - it's something I've played with from time to time.
% 
Part of the appeal of the radial basis algorithm described above is that
it is in the spline family.
                                                   Tony
---------------------------------------------------------------------
Professor A.J. Roberts     
Dept of Mathematics & Computing     E-mail: aroberts@usq.edu.au
University of Southern Queensland   Phone:  (076) 312943
Toowoomba, Queensland 4350          Fax:    (076) 312721
AUSTRALIA                           WWW: http://www.sci.usq.edu.au
                                    /pub/MC/staff/robertsa/home.html
---------------------------------------------------------------------

There is no completely satisfactory solution to the specfication of
the complex data type in high level languages because two of the
most important criteria (speed and storage convenience) are somewhat 
conflicting:
Speed is often the highest priority because of the necessity to perform
computationally intensive transforms like the FFT. For this the FORTRAN
convention of interleaving real/imaginary is fastest because the complex
data can be localized and stays in the CPU cache longer. Separating
the real and imaginary causes more extended memory references and
flushes the cache. The speed difference can be dramatic: two or three
fold on a Sparc for large arrays (64k of doubles etc).
But interleaving means always reserving space for imaginary components
even when you don't need them. Using separate, dynamic arrays for
real and imaginary is convenient (but slower) because you don't have to
allocate the imaginary if you don't need them.
How do I resolve this conflict? I use BOTH structures. Most of the time
I keep separate structures for real and imag, but copy to FORTRAN
style structures for large FFT's and back to separate again. The time
used for copying is easily made up for in the time saved by cache
locality. The rest of the time I save memory at the expense of a
slight loss in speed.
Another tip: avoid C++ objects and stick to plain old extern C structs, even
when writing in C. This will increase your portability to other languages
(Delphi etc) via DLL's (Windows) or Unix object libraries. Keep it simple.
John Day

>Another tip: avoid C++ objects and stick to plain old extern C structs, even
>when writing in C. This will increase your portability to other languages
                        ^^^ of course I meant to say "in C++"

Joel Shellman wrote:
> 
> Has there been any good methods developed for numerical solution of
> DAE's?  Is there anywhere on the net that has information about it?  I
> read a book recently and it said there wasn't a good method for this
> yet.  The book is a few years old, so I'm wondering what the current
> situation is.
Have a look at Linda Petzold's DASSL and DASPK.  As far as performance
of numerical methods go, my understanding is that it all depends on the
index of your DAE system.  For low-index problems, these codes do pretty
well.
For more info and codes, try: http://www.cs.umn.edu/~petzold/
-- 
| Phil_Calvin@unc.edu                                                  |
| Center for Multiphase Research                                       |
| Environmental Sci. and Eng.                                          |
| UNC-Chapel Hill                                                      |

Joel Shellman wrote:
> 
> Has there been any good methods developed for numerical solution of
> DAE's?  Is there anywhere on the net that has information about it?  I
> read a book recently and it said there wasn't a good method for this
> yet.  The book is a few years old, so I'm wondering what the current
> situation is.
> 
> Thanks,
> 
> -joel
Hello Joel,
Yes, there are some. Look at the ode-directory at Netlib.
One of the most general methods is the DASSL-algorithm witch
is based on the backward differetiation formulas. This
algorithm is able to solve fully implicit DAE's with the
index one. (/pub/netlib/ode/ddassl.f, see: 'The Solution of 
Initial Value Problems in Ordinary Differential-Algebraic
Equations'; Brenan, K.E.; Campbell, S.L.; Petzold, L.R.)
Some other codes are:
 o LSODI (Hindmarsh, A.C.), linear-implicit DAE, (**)
 o LIMEX (Deuflhard, P.), linear-implicit DAE, also a sparse-version
   available, see ftp://elib.zib.de/pub/elib/codelib/limexs/limexs.f
 o In the Harwell-Library you will find the Code DC03,
   also based on the BDF-Formulas, should be able to solve
   fully implicit DAE's with index 1
 o DAESOL (Eich, E.; Bock, H.G.), also based on BDF, should be 
   able to solve fully implicit DAE's with index 1, newer code and       
   algorithm than DAASL. (**)
 o RADAU5 (Hairer, E.; Wanner, G.), linear-imlicit DAE, 
   based on an implicit Runge-Kutta-Algorithm, with some 
   caution the solution of index-3-problems is posssible,
   see: http://www.ma.ic.ac.uk/~jcash/ODE_software/radau5.f
   (**) means:  don't know where to find or where to buy, respectively.
To get more Information try altavista on the names and keywords.
Hope it helps,
               Jens Bausa.
*****************************************************************
*  Dipl.-Ing. Jens Bausa                                        *
*  TU Berlin                                                    *
*  Institut fuer Energietechnik                                 *
*  Sekretariat KT 1                                             *
*  Marchstrasse 18                                              *
*  D-10587 Berlin                                               *
*  Tel.   : (..49) 30/314-26513 (-22181)                        *
*  Fax.   : (..49) 30/314-21683                                 *
*  E-Mail : bausa@olymp.fb10.tu-berlin.de                       *
*****************************************************************

In article <554t1u$i1h@winx03.informatik.uni-wuerzburg.de>, Manuel
Kessler  writes
>
>
>Joseph Sirosh (sirosh@cs.utexas.edu) wrote:
>:> HI! Is anybody aware of a source for BLAS routines optimized
>:> for the Pentium and Pentium Pro machines? Or a signal processing
>:> library optimized for these machines? Since these are superscalar
>:> processors, I've heard that proper optimization can greatly improve 
>:> the speed of typical vector operations...
Intel have a FREE CD with their Maths Kernel libraries and a signal
processing library. This contains the optimised BLAS and FFT routines.
You can get the CD direct from Intel support. I think it is called
something like performance Beta 4 CD.
-- 
Ian Parker