mm-281
Subject: Re: Constructing columns> I have indexed only
the non-zero entries of a sparse matrix> using the expression:
ij=(row-1)*matrix_size +columns. Is it possible> to reconstruct
the the row and the columns of the matrix depending> only the
index, ij. Where matrix_size is the size of a given square>
matrix. This is because I can not save the row and the columns
of each> entry in the first pass of my program.Here's how I
would think through this:ij/matrix_size = (row-1) +
columns/matrix_sizesoINT(ij/matrix_size) = row-1thus,row = 1 +
INT(ij/matrix_size)columns = ij -
(row-1)*matrix_sizehth,Rick
===
Subject: math is for wieners
onlyReal men dont need math.They will get their +3 battle axes
and hacktheir opponents to pieces before they can calculate
anything.So better improve your fighting skill because good
math skills will beof NO USE if 10 enemies are standing before
you and want to kill you !
===
Subject: Re: math is for wieners
onlyIn sci.math.num-analysis, Alberto
Panno-Peano<moltotarl@yahoo.com.mx><a5a4ac7a
.0401172131.5d697161@posting.google.com>:> Real men dont need
math.They will get their +3 battle axes and hack> their
opponents to pieces before they can calculate anything.> So
better improve your fighting skill because good math skills
will be> of NO USE if 10 enemies are standing before you and
want to kill you !Ooookaaaay.....berserker attacks Soldier #1
with +3 battle axe: 6 HP damage!Soldier #1 is
wounded.berserker attacks Soldier #1 with +3 battle axe: 4 HP
damageSoldier #1 is gravely wounded.Soldier #2 comes
up.Soldier #2 attacks berserker with M-16. *critical hit*, 8
HP damage!!berserker is gravely wounded.Medical Attendant #1
comes up.Medical Attendant attacks Solder #1 with bandage: -1
HP damage.Medical Attendant attacks Solder #1 with neck brace:
-1 HP damage.Medical Attendant attacks Solder #1 with IV
device: -1 HP damage.Medical Attendant attacks Solder #1 with
guerney strap cart: -1 HP damage.Soldier #1 is carried off
melee battle field.Soldier #3 comes up.Soldier #3 attacks
berserker with fists, intent to capture.berserker makes saving
throw.Soldier #4 comes up.berserker attacks Soldier #3 with +3
battle axe (-1 bloodslick):fumble; axe dropped.Soldier #5
comes up.Soldier #3 attacks berserker with fists, intent to
capture.Soldier #4 attacks berserker with fists, intent to
capture.berserker captured.Or perhaps:berserker attacks Tank
#1 with +3 battle axe: Active armor explosion.Active armor
attacks berserker with explosion: 12 HP damage!!berserker
dies.Or maybe:berserker attacks Plane #1 with +3 battle axe:
Plane out of range.Plane #1 launches
SeekerMissile.SeekerMissile set to berserker.SeekerMissile
attacks berserker with Explosive Warhead: 20 HP
damage.SeekerMissile attacks ground with Explosive Warhead: 16
HP damage.SeekerMissile attacks SeekerMissile with Explosive
Warhead: 18 HP damage.berserker dies.ground shows a small
crater.SeekerMissile disintegrates.Or maybe:berserker attacks
Plane #1 with +3 battle axe: Plane out of range.Plane #1
attacks berserker with Gatling Gun: *critical hit*, 9 HP
damagePlane #1 attacks berserker with Gatling Gun: *critical
hit*, 6 HP damagePlane #1 attacks berserker with Gatling Gun:
*critical hit*, 8 HP damage[... 1,497 messages
suppressed]berserker dies.Of course there's always:beserker
attacks mathematician with +3 battle axe: 5 HP
damage.mathematician is wounded.mathematician
disengages.mathematician makes saving throw.mathematician
closes door.beserker attacks door with +3 battle axe: 4 HP
damage.door shows slit.mathematician attacks cellular phone
with dialing 911: 0 HP damage.cellular phone does not make
savings throw.cellular phone is charmed!cellular phone attacks
airwaves with transmitter: 0 HP damage.airwaves does not make
savings throw.airwaves is charmed!airwaves attacks Operator
Receiver with reception: 0 HP damage.Operator Receiver does
not make savings throw.Operator Receiver is charmed!Operator
Receiver attacks Operator with message: 0 HP damage.Operator
does not make savings throw.Message sent.Operator attacks
Operator Receiver with bulletin: 0 HP damage.Operator Receiver
does not make savings throw.Operator Receiver is
charmed!Operator Receiver attacks airwaves with bulletin: 0 HP
damage.airwaves does not make savings throw.airwaves is
charmed!airwaves attacks Cruiser #1 with bulletin: 0 HP
damage.airwaves attacks Cruiser #3 with bulletin: 0 HP
damage.airwaves attacks Cruiser #8 with bulletin: 0 HP
damage.airwaves attacks Cruiser #12 with bulletin: 0 HP
damage.Cruiser #1 does not make savings throw.Cruiser #8 does
not make savings throw.Cruiser #1 is charmed!Cruiser #12 does
not make savings throw.Cruiser #8 is charmed!Cruiser #3 does
not make savings throw.Cruiser #12 is charmed!Cruiser #1 is
moving.Cruiser #3 is charmed!Cruiser #8 is moving.Cruiser #3
is moving.Cruiser #12 is moving.Cruiser #1 attacks airwaves
with siren wail: 0 HP damage.airwaves does not make savings
throw.airwaves is charmed!airwaves attacks berserker with
siren wail: 0 HP damage.No effect.Cruiser #3 attacks airwaves
with siren wail: 0 HP damage.airwaves does not make savings
throw.airwaves is charmed!airwaves attacks berserker with
siren wail: 0 HP damage.No effect.[... 10 messages
suppressed]beserker attacks door with +3 battle axe: 5 HP
damage.door shows hole.Cruiser #3 reaches
destination.berserker turns wielding +3 battle axe to attack
(dex 9);attack next melee round.Officer #6 attacks berserker
with service revolver: missed.Officer #6 attacks wall with
service revolver: 4 HP damage.wall shows small hole.Officer #9
attacks berserker with service revolver: 4 HP damage.berserker
is wounded.Cruiser #12 reaches destination.Officer #11 attacks
berserker with service revolver: 3 HP damage.berserker is
seriously wounded.Officer #18 attacks berserker with service
revolver: 3 HP damage.berserker is gravely wounded.Cruiser #8
reaches destination.Officer #15 attacks berserker with service
revolver: 3 HP damage.berserker dies.mathematician attacks
doorknob with hand: 0 HP damage.doorknob does not make savings
throw.doorknob is charmed!door is now open.mathematician
attacks Officer #6 with thanks: -1 HP damage.Officer #6 does
not make savings throw.Officer #6 is charmed!Officer #6
attacks airwaves with communications radio: 0 HP
damage.[...]-- #191, ewill3@earthlink.net -- insert random
misspent youth hereIt's still legal to go .sigless.
===
Subject:
Need direction on solving an eigenvalue problem I am trying to
solve an eigenvalue problem, and since thereare so many linear
algebra experts in this group, I am hoping that Ican get some
direction here. I have two positive integers T and k (T>k),
and a(T-k+1)x(T-k+1) symmetric Toeplitz matrix A with its
(i,j)th elementgiven by A(i,j) = max(k-abs(i-j),0)-k^2/T i.e.,
the diagonal elements are k-k^2/T, the firstsubdiagonals are
(k-1)-k^2/T, and so on until the (k-1)th subdiagonalsare
1-k^2/T, beyond that all the elements are -k^2/T (there may
not besuch elements though if T is small) . Is there a closed
form solution to the eigenvalues (and/oreigenvectors) of this
matrix? I know for band tridiagonal symmetricToeplitz matrix,
there is a closed form solution which involves thecosine
function. But is there one for such a simple matrix which
isjust a function of T and k? If there is no analytical
solution to this problem, thenis there an efficient method to
solve this problem? This is because Tcan be very large and
without an algorithm to deal with it, it isquite difficult to
solve it using eig. Many thanks for all your help.Raymond
Kankan@chass.utoronto.ca
===
Subject: Re: Data fitting problem
with chebyshev polynomials> Hi> I want to fit some pure
complex data x (only imaginary part) to some other> complex
data, y, by a polynomial with linear least squares.> y=p0
T0(x) +p1 T1(x)+...+pn Tn(x)> where Ti(x) represents the i'th
Chebyshev polynomial.> (i thought this would provide me with
better numerical conditioning,> since i get Vandermonde system
to solve).> I need real coefficients in my polynomial, so i
assume p0, ..., pn and the coefficients of Ti(x) have to be
real,> while x is complex. I calculate coefficients with
(Ax=b)> With left matrix A :> Re(T0(x)) Re(T1(x)) ....
Re(Tn(x))> Im(T0(x)) Re(T1(x)) .... Re(Tn(x))> ... for all x>
unknown coefficients are X=[p0 p1 ... pn]' > And right
columnvector vector b :> Re(y)> Imag(y)> ... for all y> Every
equation is split in real and imaginary part to make
coefficients> real.> And there's my problem. This orthogonal
technique works well if x is real> and scaled in [0,1], since
Ti(x) doesn't return very high values (only between> 0 and 1).
But if x is complex, this property is not longer valid !> e.g.
T2(x)=2x^2-1 if x=1*i then T2(x)=-3> In fact, this method gets
every sooner illconditioned than when i don't> use Chebyshev
polynomials. > Does anyone know what i'm doing wrong, or does
this technique only work> for real data??Chebyshev is good
only if the arguments are in the interval [-1,1];the function
values may be arbitrary.If your arguments are purely
imaginary, use a linear transformation tomove them into
[-1,1]. If they are scattered complex, it is best tocompute
their mean z and use as basis for expansion the powers
(x-z)^k.If they are strongly correla complex, move the mean to
zero and thelargest singular value to some real number s in
[0,1] such that the cloudof values fits inside an ellipse with
focal point at s going through 1,and use the Chebyshev
basis.Arnold Neumaier
===
Subject: Re: Number sequences,
programming, preparing for Programming Competition> I am
preparing for a student programming competition and
from[snip]> One of the problems I don't know answer is:> A
sequence of number powers of 2 starting from 1 is given,> find
the Nth digit of the sequence.> Example:> 1 2 4 9 16 25 36 49
64 81 100 121 144 ...> 7th digit is 2, 14th is 4, etc.>
P.KruminsI guess there is no other way but to compute the
squares, as many as needed.It should be useful to keep in mind
that when computing n^2 you can use theprevious square:n^2 =
(n-1)^2 +1 +2*(n-1)and (n-1)^2 you know from previous
computation. This saves a lot of workwhen n is big.Hope this
helps,Teijo.
===
Subject: Re: Number sequences, programming,
preparing for Programming Competition I am preparing for a
student programming competition and from> [snip]> One of the
problems I don't know answer is:> A sequence of number powers
of 2 starting from 1 is given,> find the Nth digit of the
sequence. Example:> 1 2 4 9 16 25 36 49 64 81 100 121 144 ...
7th digit is 2, 14th is 4, etc.> P.Krumins> I guess there is
no other way but to compute the squares, as many as needed.I
disagree.You can calculate the boundaries where the number of
digitsin a square increases. So you can write code that looks
likethis: n := n-1 k := 1 d := 1 loop: # invariants: # - n >=
0 # - we want digit number n starting from the first # digit
of k^2, counting from 0 # - k^2 is the smallest square with >=
d digits next_k := ceiling(10^[d/2]) if d*(next_k-k) <= n: n :=
n - d*(next_k-k) d := d + 1 k := next_k go round the loop again
else: m := floor(n/d) k := k + m n := n - d*m take the n'th
digit of k^2 (counting from 0) we're doneI just did a simple
implementation of this; it takes about 1/3of a second on my
machine (1GHz Athlon) to find that the 10^400'thdigit of the
sequence is 8.-- Gareth McCaughan.sig under
construc
===
Subject: Re: Number sequences, programming,
preparing for Programming Competition>> I guess there is no
other way but to compute the squares, as many as>> needed. > I
disagree.> You can calculate the boundaries where the number of
digits> in a square increases. So you can write code that looks
like> this:> n := n-1> k := 1> d := 1> loop:> # invariants:> #
- n >= 0> # - we want digit number n starting from the first>
# digit of k^2, counting from 0> # - k^2 is the smallest
square with >= d digits> next_k := ceiling(10^[d/2])> if
d*(next_k-k) <= n:> n := n - d*(next_k-k)> d := d + 1> k :=
next_k> go round the loop again> else:> m := floor(n/d)> k :=
k + m> n := n - d*m> take the n'th digit of k^2 (counting from
0)> we're done> I just did a simple implementation of this; it
takes about 1/3> of a second on my machine (1GHz Athlon) to
find that the 10^400'th> digit of the sequence is 8.Thanks!
This was very helpful.P.Krumins
===
Subject: NMR?!!can anyone
suggest some basic reference book for NMR?thanks
===
Subject:
Combination Problem I would appreciate is anybody could help
me with the following question: lets A be such that : 0 <= A
<= 4294967296 (i.e. 2^32) what is the formula giving me the
number of 4 integers combinations (N1, N2, N3, N4) (order does
not matter) 0 <= N1,N2,N3,N4 <= 4294967296 that comply to the
following requirement: N1 + N2 + N3 + N4 = Ae.g. A = 53
(53,0,0,0) (24,13,7,9) (50,2,1,0) ......Thanks for the
help,DanSubject: Re: Combination Problem
===
> I would
appreciate is anybody could help me with the following
question:> lets A be such that : 0 <= A <= 4294967296 (i.e.
2^32)> what is the formula giving me the number of 4 integers
combinations> (N1, N2, N3, N4) (order does not matter)> 0 <=
N1,N2,N3,N4 <= 4294967296> that comply to the following
requirement: N1 + N2 + N3 + N4 = A> e.g.> A = 53> (53,0,0,0)
(24,13,7,9) (50,2,1,0) ......> Thanks for the help,> DanYou
are looking for the Partition Function q:See:
http://mathworld.wolfram.com/PartitionFunctionq.htmlHugo
Pfoertner
===
Subject: Re: Combination Problem> I would
appreciate is anybody could help me with the following
question:> lets A be such that : 0 <= A <= 4294967296 (i.e.
2^32)> what is the formula giving me the number of 4 integers
combinations> (N1, N2, N3, N4) (order does not matter)> 0 <=
N1,N2,N3,N4 <= 4294967296> that comply to the following
requirement: N1 + N2 + N3 + N4 = A> e.g.> A = 53> (53,0,0,0)
(24,13,7,9) (50,2,1,0) ......> Thanks for the help,> Dan> You
are looking for the Partition Function q:> See:
http://mathworld.wolfram.com/PartitionFunctionq.html> Hugo
PfoertnerThanks for your assistance.After reviewing the
Partition Function q web page it turned out thatit was not
exactly the function I was looking for as to my
bestunderstanding it calcula the number of partitions where
all theintegers are distinct.In my case I desired also the
cases where all the integers are notdistinct.However, this
lead me to the P(n,k) function which provided thedesired
calculation with one exception : it does not take the 0
(zero)value as a possibility (i.e. (10,0,0) is not a partition
coun byP(10,3)).Now I'm left with finding the number of
partitions containing zero.Again, thanks.DanSubject: Re:
Combination Problem
===
> I would appreciate is anybody could
help me with the following question: lets A be such that : 0
<= A <= 4294967296 (i.e. 2^32) what is the formula giving me
the number of 4 integers combinations> (N1, N2, N3, N4) (order
does not matter) 0 <= N1,N2,N3,N4 <= 4294967296 that comply to
the following requirement: N1 + N2 + N3 + N4 = Ae.g.> A = 53>
(53,0,0,0) (24,13,7,9) (50,2,1,0) ......Thanks for the help,>
DanYou are looking for the Partition Function q:See:
http://mathworld.wolfram.com/PartitionFunctionq.htmlHugo
Pfoertner> Thanks for your assistance.> After reviewing the
Partition Function q web page it turned out that> it was not
exactly the function I was looking for as to my best>
understanding it calcula the number of partitions where all
the> integers are distinct.The number of partitions of n with
k or fewer addends,.....An example is given:partitions of 5
into three or fewer parts are (5),(4,1),(3,2),(3,1,1) and
(2,2,1)There is no restriction to all integers being distinct
(..1,1),(2,2,1)> In my case I desired also the cases where all
the integers are not> distinct.> However, this lead me to the
P(n,k) function which provided the> desired calculation with
one exception : it does not take the 0 (zero)> value as a
possibility (i.e. (10,0,0) is not a partition coun by>
P(10,3)).> Now I'm left with finding the number of partitions
containing zero.> Again, thanks.> DanHugo
===
Subject: Re:
Combination Problem> I would appreciate is anybody could help
me with the following question: lets A be such that : 0 <= A
<= 4294967296 (i.e. 2^32) what is the formula giving me the
number of 4 integers combinations> (N1, N2, N3, N4) (order
does not matter) 0 <= N1,N2,N3,N4 <= 4294967296 that comply to
the following requirement: N1 + N2 + N3 + N4 = Ae.g.> A = 53>
(53,0,0,0) (24,13,7,9) (50,2,1,0) ......Thanks for the help,>
DanYou are looking for the Partition Function q:See:
http://mathworld.wolfram.com/PartitionFunctionq.htmlHugo
Pfoertner> Thanks for your assistance.> After reviewing the
Partition Function q web page it turned out that> it was not
exactly the function I was looking for as to my best>
understanding it calcula the number of partitions where all
the> integers are distinct.> The number of partitions of n
with k or fewer addends,.....> An example is given:>
partitions of 5 into three or fewer parts are >
(5),(4,1),(3,2),(3,1,1) and (2,2,1)> There is no restriction
to all integers being distinct (..1,1),(2,2,1)> In my case I
desired also the cases where all the integers are not>
distinct.> However, this lead me to the P(n,k) function which
provided the> desired calculation with one exception : it does
not take the 0 (zero)> value as a possibility (i.e. (10,0,0) is
not a partition coun by> P(10,3)).> Now I'm left with finding
the number of partitions containing zero.> Again, thanks.>
Dan> HugoHugo,My apologies, you are right.I mixed things up
with the Partition Function Q.Thanks again,Dan
===
Subject: Re:
Combination Problem>> Now I'm left with finding the number of
partitions containing zero.Much easier is to take a solution
n1+n2+... = N+kand map it into (n1-1)+(n2-1)+... =
N
===
Subject: how to generate correlate matrix with specified
eigenvalues??Hi all.I am wondering if it is possible to
generate a symmetric positivematrix with 1's on the diagonals,
and with eigenvalues that I canpre-specify.If I do (in
matlab)u=rand(3,3); % 3 by 3 matrix of u(0,1) H=orth(u); %
orthonormal basis of uD=diag([.4 0 0]) the true eigenvalues
are .4 0 0x=H*D*H' % covariance matrixand then I do svd(x),
the eigenvalues of x are .4,0,0 as expec.Is there some way to
do a similar for the correlation matrix (with 1on the
diagonals) and still allow me to specify the eigenvalue?Thanks
in advance.
===
Subject: Re: how to generate correlate matrix
with specified eigenvalues?? > I am wondering if it is
possible to generate a symmetric positive> matrix with 1's on
the diagonals, and with eigenvalues that I can> pre-specify.>
If I do (in matlab)> u=rand(3,3); % 3 by 3 matrix of u(0,1) >
H=orth(u); % orthonormal basis of u> D=diag([.4 0 0]) the true
eigenvalues are .4 0 0> x=H*D*H' % covariance matrix> and then
I do svd(x), the eigenvalues of x are .4,0,0 as expec.> Is
there some way to do a similar for the correlation matrix
(with 1> on the diagonals) and still allow me to specify the
eigenvalue?Yes, by reversing the Jacobi process for finding
eigenvalues:R = ...(T3(T2(T1(D)T1')T2')T3')...,where D is the
diagonal matrix of eigenvalues, and each T is a single-plane
rotation matrix [i.e, an identity matrix,except for the
elements in positions (i,i), (i,j), (j,i), (j,j),which have
values cos a, sin a, -sin a, cos a)].The angle of rotation is
chosen so that the values inpositions (i,i) & (j,j) of T()T'
will equal one another.At convergence, the accumula product
...(T3(T2(T1))) will be thematrix of eigenvectors. R will not
be unique, but will depend on thesequence in which the (i,j)
pairs are rota.
===
Subject: Re: how to generate correlate
matrix with specified eigenvalues??> I am wondering if it is
possible to generate a symmetric positive> matrix with 1's on
the diagonals, and with eigenvalues that I can> pre-specify.>
If I do (in matlab)> u=rand(3,3); % 3 by 3 matrix of u(0,1) >
H=orth(u); % orthonormal basis of u> D=diag([.4 0 0]) the true
eigenvalues are .4 0 0> x=H*D*H' % covariance matrix> and then
I do svd(x), the eigenvalues of x are .4,0,0 as expec.> Is
there some way to do a similar for the correlation matrix
(with 1> on the diagonals) and still allow me to specify the
eigenvalue?> Yes, by reversing the Jacobi process for finding
eigenvalues:> R = ...(T3(T2(T1(D)T1')T2')T3')...,> where D is
the diagonal matrix of eigenvalues, > and each T is a
single-plane rotation matrix [i.e, an identity matrix,> except
for the elements in positions (i,i), (i,j), (j,i), (j,j),>
which have values cos a, sin a, -sin a, cos a)].> The angle of
rotation is chosen so that the values in> positions (i,i) &
(j,j) of T()T' will equal one another.So the diagonal entries
will equal one another -- but it is not possible, in general,
to generate a symmetric positive matrix with 1's on the
diagonal having specified eigenvalues, as reques. For example,
consider the most general symmetric 2 x 2 matrix with 1's on
the diagonal: 1 a a 1This matix has eigenvalues l1 = 1-a and
l2 = a+1 (and is positive definite for -1 < a < 1). Hence,
only if the specified eigenvalues satisfy l1 + l2 = 2 can a
matrix of the reques form be found.Paul-- Paul Abbott Phone:
+61 8 9380 2734School of Physics, M013 Fax: +61 8 9380 1014The
University of Western Australia (CRICOS Provider No 00126G) 35
Stirling HighwayCrawley WA 6009 mailto:paul@physics.uwa.edu.au
AUSTRALIA http://physics.uwa.edu.au/~paul
===
Subject: Re: how to
generate correlate matrix with specified eigenvalues??> [...]>
So the diagonal entries will equal one another -- but it is
not > possible, in general, to generate a symmetric positive
matrix with 1's > on the diagonal having specified
eigenvalues, as reques. For example, > consider the most
general symmetric 2 x 2 matrix with 1's on the > diagonal:> 1
a> a 1> This matix has eigenvalues l1 = 1-a and l2 = a+1 (and
is positive > definite for -1 < a < 1). Hence, only if the
specified eigenvalues > satisfy l1 + l2 = 2 can a matrix of
the reques form be found.Yes, the sum of the specified
eigenvalues must equal the order of thematrix, or they can not
possibly be the eigenvalues of a correlation matrix. I took
that as given.
===
Subject: Re: how to generate correlate matrix
with specified eigenvalues??> I am wondering if it is
possible to generate a symmetric positive> matrix with 1's on
the diagonals, and with eigenvalues that I can> pre-specify.>
If I do (in matlab)> u=rand(3,3); % 3 by 3 matrix of u(0,1) >
H=orth(u); % orthonormal basis of u> D=diag([.4 0 0]) the true
eigenvalues are .4 0 0> x=H*D*H' % covariance matrix> and then
I do svd(x), the eigenvalues of x are .4,0,0 as expec.> Is
there some way to do a similar for the correlation matrix
(with 1> on the diagonals) and still allow me to specify the
eigenvalue?> Yes, by reversing the Jacobi process for finding
eigenvalues:> R = ...(T3(T2(T1(D)T1')T2')T3')...,> where D is
the diagonal matrix of eigenvalues, > and each T is a
single-plane rotation matrix [i.e, an identity matrix,> except
for the elements in positions (i,i), (i,j), (j,i), (j,j),>
which have values cos a, sin a, -sin a, cos a)].> The angle of
rotation is chosen so that the values in> positions (i,i) &
(j,j) of T()T' will equal one another.> At convergence, the
accumula product ...(T3(T2(T1))) will be the> matrix of
eigenvectors. R will not be unique, but will depend on the>
sequence in which the (i,j) pairs are rota.Thanks for the tip.
Would you have a reference or an example I canunderstand better
how to do this? How would I find the angle ofroation, and what
is 'a' in cos(a)??Thanks again.
===
Subject: Re: how to generate
correlate matrix with specified eigenvalues?? S. Eng schrieb::>
Thanks for the tip. Would you have a reference or an example I
can> understand better how to do this? How would I find the
angle of> roation, and what is 'a' in cos(a)??> Thanks
again.If you have a correlation-matrix R, then you can bring
itto a similar diagonal-matrix D by rotations, say T.How to
get a proper rotation matrix may be a second step fornow. If T
is a rotation.matrix, then it is orthogonal and T * T' = INow
if you find a certain T, given C, which rotates itto a similar
diagonalmatrix D you write this in a formula D = T*R*T'T is
then the rotation, which finds the eigenvalues ...it's just
the eigenvectors of R.Normally you *HAVE* an R and find, for
instance bycholesky-decomposition of R and seeking the
PCA-positionof the loadingsmatrix a rotation-matrix T
(implicitely)by the software. Normally the software doesn't
give accessto this rotation-matrix, so this is not well known,
usually.But now your case is opposite: you *have* a set of
eigenvaluesand want to create a correlation-matrix. Just
create any arbitraryrotation-matrix, and the appropriate
correlation-matrix willbe genera.Now, how to create any
rotation-matrix? It's just to multiplytwo simpler
rotation-matrices. For instance, you guess an angle ofsay,
0.2*pi, let say applied to the pair of the first two
componentsin D. Then say, c is the cosine(0.2*pi) and s is the
sine(0.2*pi).All other components stay untouched. If you have 4
dimensions, thenyour T1 is then T1 =
((c,s,0,0)(-s,c,0,0),(0,0,1,0),(0,0,0,1))First rotation is
then TMP = T1' * D * T1Then you select another pair of
components and another angleto rotate. That's just the
opposite way, that Jacobi invenhis algorithm to find the
principal components position, givenR.In the end you have
(using all possible pairs) R = T6'*T5'*T4'*T3'*T2'*T1' * D *
T1*T2*T3*T4*T5*T6and you could have crea a random
rotationmatrix as a wholeby simply multiplying T =
T1*T2*t3*T4*T5*T6 and R = T' * D * TT is then the eigenvector
of R, btw, and D^(1/2)*T is just thePCA-loadingsmatrix P of R,
with the properties: P *P' = R P'*P = Dwhich will be found, if
you put R into your PCA-analysis program.Gottfried
Helms
===
Subject: Re: how to generate correlate matrix with
specified eigenvalues??> [...] Would you have a reference or
an example I can> understand better how to do this? How would
I find the angle of> roation, and what is 'a' in cos(a)??'a'
is the angle of rotation. Since you didn't recognize this,
Iconclude that you also are not familiar with the Jacobi
method forgetting eigenvalues. You should study it -- it
should be covered inany good text on numerical methods for
matrices (e.g., Golub & VanLoan's_Matrix Computations_) -- and
then write your own Jacobi routine.Once you have it working,
you should have no trouble modifying itto run backwards,
starting from a diagonal matrix of eigenvaluesand finishing
when all the diagonal elements of R are equal.To rotate in the
(i,j) plane, let p,r,r,q denote the valuesin positions
(i,i),(i,j),(j,i),(j,j) of the current R.To set the new values
in positions (i,i) & (j,j) to (p+q)/2,set tan(2a) =
(q-p)sgn(r)/|2r|, with 2a in [-Pi/2,Pi/2];adapt the half-angle
reduction formula tan(a) = sin(2a)/(1+cos(2a))to get sin(a) &
cos(a), with a in [-Pi/4,Pi/4].
===
Subject: Re: PI> Hi> In a
book I have found the following formulation[...]> A capital
'pi' means product, just like a capital 'sigma' means
summation.> OUP
===
==Two rational approximations of ,,pi
.
===
====================================Suppose that c is in
the set {-1/4,1/4} and define n^2+(n-c)(2n+3-4c)
n^2(n-2c)^2a_n(c)= ------------------ ,
b_n(c)=---------------------------- , n=0,1,... .
2(n-c)(n-c+1) 4(2n-1-2c)(2n+1-2c)(n-c)^2 Consider the
sequences (A_n(c))_{n>=0} and (B_n(c))_{n>=0} where
A_{n+1}(a)=a_n(c)*A_{n}(c) -b_n(c)*A_{n-1}(c)(1) , n=1,2,....,
B_{n+1}(a)=a_n(c)*B_{n}(c) -b_n(c)*B_{n-1}(c)with 4(17-8c)
A_0(c)=2(1-4c) , A_1(c)=-------- 15 2(11-4c) B_0(c)=1 ,
B_1(c)=--------- . 15 Define A_n(1/4) A_n(-1/4)(2) q_n=
--------- , Q_n= ----------- , n=0,1,... . B_n(1/4) B_n(-1/4)
Then(3) 3=q_1<...<q_n<q_{n+1}<...< pi < ...< Q_{k+1}< Q_k<...<
Q_1=19/6 .Moreover, there exists C_1>0 , C_2>0 such that for
all natural numbers n C_2*32^{-n} 0 < pi- q_n < C_1*32^{-n} ,
o< Q_n - pi < ------------ . sqrt(n)Because (1) si simple,
numbers q_n and Q_n from (2) are easy to be obtained. Further
use (3). I appreciate that the last inequalitiesare not the
best.
===
====================Subject: Re: PI> HiIn a book I
have found the following formulation[...]> A capital 'pi'
means product, just like a capital 'sigma' means summation.>
OUP> 
===
==> Two rational approximations of ,,pi
.[............]> 
===
====================================Some
records :1949 = John von Neumann using computer ENIAC has
found the first 2037decimals of pi.1973= Guillound and Bouyer,
using CDC-7600 together with someMachin-like formulas, present
the first 1,000,000 decimals ofpi1983= Y. Kanada, Y. Tamura,
S.Yoshiro and Y. Ushiro find the first 10,013,395 decimals.
1989= David Chudnovski and Gregory Chudnovski have put in
evidence 1,000,000,000 of pi 1988= D.H.Bailey ,,The
Computation of Pi to 29,360,000 Decimal DigitesUsing Borwein's
Quartically Convergent Algorithm, Mathematics of Computation ,
(1988) 283-296. 1995= Simon Plouffe, Peter B. Borwein and D.H.
Bailey have announcedthat the 1010-hexa digit of pi is 9. 2000
(March 25)= Colin Martin - Sydney, Australia found
1,610,612,736 (1.5-gig) digits. using WindowsNT4.0 found
17,100,000,000 digits(15.93-gigs)Some incomplete references ,
sorry for misprints, are : { Please see [5] , [6-B], [29] }
[1] D. H. Bailey}, ,,The Computation of pi to 29.360.000
DecimalDigits Using Borwein s Quartically Convergent
Algorithm, Mathematical Computations, nr. 50, pg
283-296,1988[2] D. H. Bailey, P. Borwein and Simon Plouffe},
,,On the RapidComputation of Various Polylogarithmic
Constants,1996,see:
http://www.cecm.sfu.ca/~pborwein/PAPERS/PI23.ps[3] David H.
Bailey, Johnatan M. Borwein and Simon Plouffe, ,,TheQuest for
Pi , The Mathematical Intelligencer, June 1996[4] David H.
Bailey, ,,A Compendium of BBP - Type Formulas forMathematical
Constants,( 2000)[5] Fabrice Bellard, ,,Computation of the
n-th digit of ,,pi in anybase inbig{O}(n^{2}, (1997) (try to
find on Internet)[6] Jonathan M Borwein and Peter Borwein,
,,Means, Iterations and Experimentaly Induced Identities,
Simon Fraser University[6-B] Jonathan M. Borwein , ,,The
Experimental Mathematician: AComputational Guide to the
Mathematical Unknown , (2002) (nice !),see:
www.cec,.sfu.ca/~jborwein/talks.html [7] Kurt Van den Brenden
(Branden ?) , ,,Fun with Pi[8] Steven Finch}, ,,Archimedes
Constant, see: http://www.mathsoft.com/constants[9] Steven
Finch, ,,The constant ,, pi ...The Clasic period , 11Noiembrie
1999, see: http://www.mathsoft.com/piclas~1.html [10] Steven
Finch, ,,Series formulas for ,, pi , 5 Martie 2000 ,see
http://www.mathsoft.com/piseries.html[11] Steven Finch, ,,The
Miraculous Bailey - Borwein - Plouffe PiAlgorithm, see:
http:// www.mathsoft.com/BBPalgorithm.html[12] Steven Finch,
,, The constant pi ...The geometric period,(2000), see:
http://www.mathsoft.com/PIGEOM~1.html[13] Steven Finch, ,, The
constant $pi$ ... 4 .... Pi and theArithmetic - Geometric Mean,
(2000), see : http://www.mathsoft.com/PIAGM~1.html[14] Steven
Finch, ,,N-thdigit computation , 15 Feb 2000, see:
www.mathsoft.com/NTHDIGIT.htm[15] Steven Finch , ,,The
Miraculous Bailey - Borwein - Plouffe PiAlgorithm, 1 Oct 1995,
see: http://www.mathsoft.com/TheBBPAlgorithm.htm [16] Steven
Finch, ,,Numerical Algorithms , 17 Aug 2000, see:
http://numbers.computation.free.fr/constants/Pi/
iterativePi.html[17] R. Fueter, ,,Vorlesungen uber die
singularen moduln und diekomplexe multiplikation der
elliptichen Funktionen, 1926,Berlin.[18] B. Gourevitch ,
,,John Machin, see:
http://www.multimania.com/bgourewitch/mathematiciens/machin/
machi.html[19] B. Gourevitch, ,,Formule dite de BBP, Bailey -
Borwein -Plouffe,
see:http://multimania.com/bgourevitch/mathematiciens/plouffe/
plouffe.html[20] M. Hirschhorn, ,,A new formula for ,, pi , 5
Dec 1997 [21] Ron Knott, ,,Proof of an Arctan formula ,(1998),
see:
http://www.mcs.surrey.ac.uk/Personal/R.Knott/Fibonacci/
arctanproof1.html[22] Ron Knott, ,,Pi and the Fibonacci
numbers},( 1999) , see
http://www.mcs.surrey.ac.uk/Personal/R.Knott/Fibonacci/
fibpi.html[23] Lazarus Mudehwe, ,,The story of Pi ... [24] J.
J. OConnor and E.F. Robertson, ,,A history of ,, pi ,(2001),
see
http://wwwww-history.mcs.st_andrews.ac.uk/history/HistTopics/
Pi_through_the_ages.html[25] Simon Plouffe, ,,Table of current
records for the computaation ofconstants,(1999), see:
http://www.lacim.uqam.ca/pi/records.html[26] Simon Plouffe,
,,1 Billion Digits of Pi,see:
http://www.lacim.uqam.ca/piDATA/PI.html[27] Simon Plouffe,
,,On the computation of n-th decimal digit ofvarious
transcendental numbers,(1996), see:[28] S. Rabinowitz}, ,,A
Spigot-Algorithm for ,, pi, Abstract of the American
Mathematical Society, (1991), vol. 12,pg. 30[29] Francis Su,
,,Finding the N-th digit of Pi,....[30] Eric W. Weisstein,
,,Pi Formulaas, see:
www.mathworld.wolfram.com/PiFormulas.html[31] Eric W.
Weisstein, ,,Machin-like formulas, see:
www.mathworld.wolfram.com/Machin-LikeFormulas.html[32] Eric W.
Weisstein, ,,Pi, see: http://mathworld.wolfram.com/Pi.html[33]
Eric W. Weisstein , ,,Karatsuba Multiplication, see:
http://mathworld.wolfram.com/KaratsubaMultiplication.html[34]
Eric W. Weisstein, ,,Pi Digits, see:
http://mathworld.wolfram.com/PiDigits.html [35] ?? ,,On the
computation of the n-th decimal digit of varioustranscendental
numbesee:[36] ?? ,,Machins Merit,
http://mathpages.com/home/kmath373.html
===
====================
==Subject: a problem rela to plotting a CDF expression in
Matlab I have a question that has to do with a problem
encountered while using Matlab to verify a cumulative
distribution function (CDF) that I worked out. Hope someone on
this forum may give me some help. Anyway, here it goes. I have
derived a certain CDF expression for a random variable t, and
it is obtained from a double-integral expression. t is a
function of two variables, v and phi, and it should be able to
take on all its validity. What I got was a curve that looks
like the CDF when t is greater than a certain value (t_a). But
when t<t_a, the expression breaks down (symptom: a spike that
shoots up above 1, then comes down). I did some investigation
and believe I found the reason. In the double-integral
expresssion, the outer integral is w.r.t. phi, and the inner
integral is with respect to v. The upper and lowerbounds for
the inner integral (i.e., for v) are as follows (when
double-integrating, the range of phi is split into two):for
3*pi/2 < phi < 7*pi/4: 2Rcos(phi)/t < v < -2b/sin(phi) [1]for
7*pi/4 < phi < 2*pi: 2Rcos(phi)/t < v < 2b/cos(phi) [2]In the
above, R and b are two positive constants, and I assign 250
and 80 to them, respectively. The cause of the error is that
given R=250, b=80, when t is very small(e.g., 0.1), the lower
bound is greater than the upper bound, thus causing the error.
By using Eq. 1 and 2 above, I can also determine at what range
of t the derived CDF expression breaks down. For Eq. 1, it
breaks down when t<-(R/b)*sin(phi)*cos(phi), and for Eq. 2, it
breaks down for t<(R/b)*cos(phi)^2. So I know what the cause of
the error is, but don't know how to fix it.Since I do not wish
to place any constraints on R and b other than that they both
are positive real numbeI still would like to obtain just one
expression that takes care of the entire range of t (i.e.,
t>0) instead of having to treat the cases below and above t_a
separately. Has anybody ever encountered this kind of problem
in his/her career? I certainly would appreciate some pointer
on this. Hope to hear from you soon. Thanks.-EdSubject: Re:
SymbMath.com: web-based computer algebra system
===
> Would you
please show these results?Ah Dr Huang, thanks for addressing
my concerns.Hey I tried some of the functions and they gave
results inconsistent> with matlab. Do you think I should tell
them to fix their program?> www.SymbMath.com
===
Subject:
Explicit prime counting formula, complexityThe way I calculate
explicit prime formulas is straightforward, as Istar
withdS(N,3) = floor(N-4)/6. with even Nand itera.SoN/2 -
floor((N-4)/6) - floor((N-16)/10) + floor((N-16)/30)
-floor((N-36)/14) + floor((N-22)/42) + floor((N-106)/70)
-floor((N-106)/210) + 2.is the result of two iterations
covering 5 and 7.I'm curious now about whether or not there
are smaller expressionsthat do what it does out to infinity as
I sat down and figured out howterms get added.In general, with
x=sqrt(N) there are approximately 2^{x/ln x}/ln
Nterms.
===
Subject: Re: Explicit prime counting formula,
complexity> The way I calculate explicit prime formulas is
straightforward, as I> star with> dS(N,3) = floor(N-4)/6. with
even N> and itera.> So> N/2 - floor((N-4)/6) - floor((N-16)/10)
+ floor((N-16)/30) -> floor((N-36)/14) + floor((N-22)/42) +
floor((N-106)/70) -> floor((N-106)/210) + 2.> is the result of
two iterations covering 5 and 7.> I'm curious now about whether
or not there are smaller expressions> that do what it does out
to infinity as I sat down and figured out how> terms get
added.All these problems that are so hard on your poor little
brain are nicely solved on my webpage, together with (shall I
mention it? ) an implementation of a prime counting algorithm
that is about 1000 times faster than anything you've ever
crea...> In general, with x=sqrt(N) there are approximately
2^{x/ln x}/ln N> terms.Wrong. I mean massively, completely
wrong.
===
Subject: Re: Explicit prime counting formula,
complexity> The way I calculate explicit prime formulas is
straightforward, as I> star with> dS(N,3) = floor(N-4)/6. with
even N> and itera.> So> N/2 - floor((N-4)/6) - floor((N-16)/10)
+ floor((N-16)/30) -> floor((N-36)/14) + floor((N-22)/42) +
floor((N-106)/70) -> floor((N-106)/210) + 2.> is the result of
two iterations covering 5 and 7.> I'm curious now about whether
or not there are smaller expressions> that do what it does out
to infinity as I sat down and figured out how> terms get
added.> In general, with x=sqrt(N) there are approximately
2^{x/ln x}/ln N> terms.Should be that with x=sqrt(N), there
are approximately 2^{x/ln x - 1)terms, which I think is a lot
forcing me to rethink my assertionsbefore that the formula is
the smallest possible.
===
Subject: Re: Explicit prime counting
formula, complexity> The way I calculate explicit prime
formulas is straightforward, as I> star withdS(N,3) =
floor(N-4)/6. with even Nand itera.SoN/2 - floor((N-4)/6) -
floor((N-16)/10) + floor((N-16)/30) -> floor((N-36)/14) +
floor((N-22)/42) + floor((N-106)/70) -> floor((N-106)/210) +
2.is the result of two iterations covering 5 and 7.I'm curious
now about whether or not there are smaller expressions> that do
what it does out to infinity as I sat down and figured out how>
terms get added.In general, with x=sqrt(N) there are
approximately 2^{x/ln x}/ln N> terms.> Should be that with
x=sqrt(N), there are approximately 2^{x/ln x - 1)> terms,
which I think is a lot forcing me to rethink my assertions>
before that the formula is the smallest possible.>
Unfortunately, you have previously repor that the 'sqrt' has
an *inherent* ambiguity which isinescapable and which requires
the management of negative values in its use. You have
certainly not done soabove, and negative values of 'x' will
not work with your formulation. Hence the above fails to meet
yourown standards of acceptability. The demands of honesty and
integrity demand that you withdraw your primecounting
method.The Harris Posting Model:1) Ready!2) Position
cross-hairs squarely on foot!3) Fire!--There are two things
you must never attempt to prove: the unprovable -- and the
obvious.----http://www.crbond.com
===
Subject: Re: Explicit
prime counting formula, complexity>>The way I calculate
explicit prime formulas is straightforward, as I>>star
with>>dS(N,3) = floor(N-4)/6. with even N>>and
itera.>>So>>N/2 - floor((N-4)/6) - floor((N-16)/10) +
floor((N-16)/30) ->>floor((N-36)/14) + floor((N-22)/42) +
floor((N-106)/70) ->>floor((N-106)/210) + 2.>>is the result
of two iterations covering 5 and 7.>>I'm curious now about
whether or not there are smaller expressions>>that do what it
does out to infinity as I sat down and figured out how>>terms
get added.>>In general, with x=sqrt(N) there are
approximately 2^{x/ln x}/ln N>>terms.> Should be that with
x=sqrt(N), there are approximately 2^{x/ln x - 1)> terms,
which I think is a lot forcing me to rethink my assertions>
before that the formula is the smallest possible.Take your
time.Gib
===
Subject: Re: Explicit prime counting formula,
complexity---snip---> I'm curious now about whether or not
there are smaller expressions> that do what it does out to
infinity as I sat down and figured out how> terms get added.>
In general, with x=sqrt(N) there are approximately 2^{x/ln
x}/ln N> terms.> I thought I had an answer for you. Alas, I do
not. My primary mathematics resource, a sci.math post by a
highly regardedmember of the sci.math community, is called a
Quick Math Guide. Surprisingly, it offered no guidance to your
particular question. Possibly you can e-mail the author in
private (and keep it private).Less surprisingly, I am being
sarcastic here. The quick math guidewon't help anyone,
ever.Have a nice day,Jay
===
Subject: Re: Explicit prime
counting formula, complexityIsn't it remarkable how populart
prime counting has become with JSH since he was show to be so
spectacularly wrong about properties of the ring of algrabraic
integers?
===
Subject: Parallel Algorithms Book?I'm giving a
senior-level course next quarter in parallel computing.We'll
be using a Beowolf cluster. Any suggestions for a book
thatfocuses on numerical algorithms for parallel machines but
stillassumes no background in parallel computing (and hence
addresses basicideas of architecture, etc.)? The students will
have a fair backgroundin the basics of computing, incl.
architecture and OS ideas.I can find books focusing on
architectures, or on computer scienceproblems, or even fluid
problems, but I want an Intro. to Num.Analysis--The Parallel
Version book.
===
Subject: Re: Parallel Algorithms Book?Jeffery
J. Leader (This Address is Rarely Read)
<JeffLeader@HotMail.comI can find books focusing on
architectures, or on computer science> problems, or even fluid
problems, but I want an Intro. to Num.> Analysis--The Parallel
Version book.V.-- email: lastname at cs utk eduhomepage: cs
utk edu tilde lastname
===
Subject: Too many iterations in
tqli.cX-AUTHid: minutsilHi allI'm trying to compute
eigenvalues with tqli.c from recipes in C.I'm pretty sure I'm
doing it correctly. I get similar results withmatlab.
Occasionally though, for some matrices I get the Too
manyiterations error. Does anybody know what might lead to
that? Illconditioned matrix maybe? Any workaround?
Thanks!
===
Subject: Re: Too many iterations in tqli.cI would
always suggest using LAPACK (or CLAPACK if you prefer C)
fromwww.netlib.org as opposed to anything from Numerical
Recipes. Matlab matrixcomputations (starting in Matlab 6.0)
are performed using LAPACK. AndLAPACK is free!John> Hi all>
I'm trying to compute eigenvalues with tqli.c from recipes in
C.> I'm pretty sure I'm doing it correctly. I get similar
results with> matlab. Occasionally though, for some matrices I
get the Too many> iterations error. Does anybody know what
might lead to that? Ill> conditioned matrix maybe? Any
workaround?> Thanks!
===
Subject: Re: Too many iterations in
tqli.cX-AUTHid: minutsil> I would always suggest using LAPACK
(or CLAPACK if you prefer C) from> www.netlib.org as opposed
to anything from Numerical Recipes. Matlab matrix>
computations (starting in Matlab 6.0) are performed using
LAPACK. And> LAPACK is free!> JohnI know about lapack and
clapack, and I know they are high performancelibraries, but
they are just so cumbersome. That's a pain, for people who are
not numerical-analysts, and only use these librariesas
black-boxes.I actually write code in c++ and I know about
lapack++, and I just saw it was superseeded by TNT. That's
nice, I thought, but then I see they have their own
implementation of 1D array, complex, etc. Why do that, when
there's already vector, valarray, and complex inC++? I do have
the freedom to choose what routines I use, but is therea C++
numerical library which combines the performance of lapack
orNR, with the elegance and convenience of C++? Also, one
which doesnot re-invent the wheel and which is
free.
===
Subject: Re: Too many iterations in tqli.c >I would
always suggest using LAPACK (or CLAPACK if you prefer C) from
>www.netlib.org as opposed to anything from Numerical Recipes.
Matlab matrix >computations (starting in Matlab 6.0) are
performed using LAPACK. And >LAPACK is free! >John >> Hi
all > I'm trying to compute eigenvalues with tqli.c from
recipes in C. >> I'm pretty sure I'm doing it correctly. I get
similar results with >> matlab. Occasionally though, for some
matrices I get the Too many >> iterations error. Does anybody
know what might lead to that? Ill >> conditioned matrix maybe?
Any workaround? > Thanks! >>we had his already repea
time. the code in NRC is in principle o.k. but theyuse an
overly stringent termination criterion which assumes a small
relative errorin the elements to be reduced to zero. check
this with the error criteria inLAPACK (or the original
Wilkinson_Reinsch volume). adding a small absolute error
tolerance to the termination bound removes the
problemhthpeter
===
Subject: Re: Too many iterations in
tqli.cX-AUTHid: minutsil>I would always suggest using LAPACK
(or CLAPACK if you prefer C) from>www.netlib.org as opposed to
anything from Numerical Recipes. Matlab matrix>computations
(starting in Matlab 6.0) are performed using LAPACK.
And>LAPACK is free!>>John> Hi allI'm trying to compute
eigenvalues with tqli.c from recipes in C.>> I'm pretty sure
I'm doing it correctly. I get similar results with>> matlab.
Occasionally though, for some matrices I get the Too many>>
iterations error. Does anybody know what might lead to that?
Ill>> conditioned matrix maybe? Any workaround?Thanks!> we
had his already repea time. the code in NRC is in principle
o.k. but they> use an overly stringent termination criterion
which assumes a small relative error> in the elements to be
reduced to zero. check this with the error criteria in> LAPACK
(or the original Wilkinson_Reinsch volume). adding a small
absolute error > tolerance to the termination bound removes
the problem> hth> peterExactly where was this discussed? I'd
like to check it out.
===
Subject: Re: Too many iterations in
tqli.c> I would always suggest using LAPACK (or CLAPACK if you
prefer C) from> www.netlib.org as opposed to anything from
Numerical Recipes. Matlab matrix> computations (starting in
Matlab 6.0) are performed using LAPACK. And> LAPACK is free!>
JohnSome of us (particularly those who undertake to deliver
working codes tonon-experts) have a need for eigensystem and
rela routines thatoperate on matrices contained in customized
data structures and whichmust be seamlessly integra with the
rest of the code (i.e. no3rd-party DLLs). In such cases,
adapting the original source code isthe safest route. For
eigensystems, the source code can be directlytraced to the
Handbook on Linear Algebra with it's spaghetti-coded Algol68,
through Eispack (the Fortran translation). I have never
noticedany significant changes in any of these codes, whether
sold by IMSLor freely distribu by netlib. (To be fair to Drs.
Wilkinson andReinsch, they were working before the development
of structuredprogramming.)Press, et al., have at least
performed the public service of unrollingall those GOTO direc
loops into standardized control structures.I only wish they
had done so for the whole set of Handbook routines.- Bill
Frensley
===
Subject: Re: Too many iterations in tqli.c> Hi all>
I'm trying to compute eigenvalues with tqli.c from recipes in
C.> I'm pretty sure I'm doing it correctly. I get similar
results with> matlab. Occasionally though, for some matrices I
get the Too many> iterations error. Does anybody know what
might lead to that? Ill> conditioned matrix maybe? Any
workaround? If it gives the right answer when it converges,
but sometimes fails toconverge, it sounds like there might be
an error in the calculation ofthe implicit shift. If you typed
in the code yourself, carefully checkthe lines just below the
if statement that generates that error, andbe keenly aware of
the difference between the lower case letter L andthe numeral
one in the array indices.Also, as a test, you might want to
print out the total number of iterations for the cases that do
converge. Usually implicit QL onlytakes 3 or 4 iterations for
each eigenvalue.- Bill Frensley
===
Subject: Re: Too many
iterations in tqli.cX-AUTHid: minutsil>> Hi all>> I'm trying
to compute eigenvalues with tqli.c from recipes in C.>> I'm
pretty sure I'm doing it correctly. I get similar results
with>> matlab. Occasionally though, for some matrices I get
the Too many>> iterations error. Does anybody know what might
lead to that? Ill>> conditioned matrix maybe? Any workaround?
> If it gives the right answer when it converges, but
sometimes fails to> converge, it sounds like there might be an
error in the calculation of> the implicit shift. If you typed
in the code yourself, carefully check> the lines just below
the if statement that generates that error, and> be keenly
aware of the difference between the lower case letter L and>
the numeral one in the array indices.> Also, as a test, you
might want to print out the total number of > iterations for
the cases that do converge. Usually implicit QL only> takes 3
or 4 iterations for each eigenvalue.> - Bill FrensleyI prin
out the number of iterations for each index l of the outer
loop. I see that iter is indeed small, except for the
firstone, which is way over the NR bound of 30. Also, when
iter >= 30I don't exit, and the algorithm does terminate, and
after sortingthe eigenvalues, I get a difference of about
10^(-5) compared withmatlab. So I'm very temp to let it
iterate as much as it wants,but, obviously, I have no
guarantee that it will converge nor thatit produces the right
results. l: 0 ITER=185 l: 1 ITER=0 l: 2 ITER=9 l: 3 ITER=6 l:
4 ITER=5 l: 5 ITER=3 l: 6 ITER=4 l: 7 ITER=5 l: 8 ITER=2 l: 9
ITER=5 l: 10 ITER=4 l: 11 ITER=3 l: 12 ITER=5 l: 13 ITER=3 l:
14 ITER=5 l: 15 ITER=4 l: 16 ITER=3 l: 17 ITER=4 l: 18 ITER=3
l: 19 ITER=4 l: 20 ITER=4 l: 21 ITER=4 l: 22 ITER=3 l: 23
ITER=3 l: 24 ITER=3 l: 25 ITER=4 l: 26 ITER=2 l: 27 ITER=3 l:
28 ITER=4 l: 29 ITER=4 l: 30 ITER=3 l: 31 ITER=4 l: 32 ITER=4
l: 33 ITER=4 l: 34 ITER=2 l: 35 ITER=3 l: 36 ITER=3 l: 37
ITER=3 l: 38 ITER=3 l: 39 ITER=3 l: 40 ITER=4 l: 41 ITER=4 l:
42 ITER=4 l: 43 ITER=3 l: 44 ITER=2 l: 45 ITER=4 l: 46 ITER=2
l: 47 ITER=3 l: 48 ITER=4 l: 49 ITER=2 l: 50 ITER=3 l: 51
ITER=3 l: 52 ITER=5 l: 53 ITER=2 l: 54 ITER=3 l: 55 ITER=4 l:
56 ITER=4 l: 57 ITER=2 l: 58 ITER=4 l: 59 ITER=2 l: 60 ITER=3
l: 61 ITER=5 l: 62 ITER=3 l: 63 ITER=4 l: 64 ITER=2 l: 65
ITER=3 l: 66 ITER=3 l: 67 ITER=3 l: 68 ITER=3 l: 69 ITER=4 l:
70 ITER=2 l: 71 ITER=3 l: 72 ITER=3 l: 73 ITER=3 l: 74 ITER=1
l: 75 ITER=1 l: 76 ITER=2 l: 77 ITER=3 l: 78 ITER=1 l: 79
ITER=1 l: 80 ITER=0 l: 81 ITER=2 l: 82 ITER=1 l: 83 ITER=0 l:
84 ITER=2 l: 85 ITER=1 l: 86 ITER=2 l: 87 ITER=2 l: 88 ITER=1
l: 89 ITER=0 l: 90 ITER=0 l: 91 ITER=0 l: 92 ITER=3 l: 93
ITER=3 l: 94 ITER=1 l: 95 ITER=2 l: 96 ITER=2 l: 97 ITER=2 l:
98 ITER=2 l: 99 ITER=1 l: 100 ITER=2 l: 101 ITER=1 l: 102
ITER=3 l: 103 ITER=2 l: 104 ITER=1 l: 105 ITER=1 l: 106 ITER=2
l: 107 ITER=1 l: 108 ITER=1 l: 109 ITER=0 l: 110 ITER=2 l: 111
ITER=2 l: 112 ITER=2 l: 113 ITER=4 l: 114 ITER=4 l: 115 ITER=1
l: 116 ITER=1 l: 117 ITER=1 l: 118 ITER=2 l: 119 ITER=3 l: 120
ITER=2 l: 121 ITER=1 l: 122 ITER=2 l: 123 ITER=1 l: 124 ITER=0
l: 125 ITER=2 l: 126 ITER=2 l: 127 ITER=2 l: 128 ITER=4 l: 129
ITER=0 l: 130 ITER=2 l: 131 ITER=3 l: 132 ITER=2 l: 133 ITER=2
l: 134 ITER=2 l: 135 ITER=2 l: 136 ITER=1 l: 137 ITER=0 l: 138
ITER=2 l: 139 ITER=2 l: 140 ITER=1 l: 141 ITER=1 l: 142 ITER=2
l: 143 ITER=2 l: 144 ITER=2 l: 145 ITER=1 l: 146 ITER=2 l: 147
ITER=1 l: 148 ITER=1 l: 149 ITER=0 l: 150 ITER=6 l: 151 ITER=0
l: 152 ITER=4 l: 153 ITER=4 l: 154 ITER=5 l: 155 ITER=1 l: 156
ITER=2 l: 157 ITER=1 l: 158 ITER=5 l: 159 ITER=1 l: 160 ITER=0
l: 161 ITER=2 l: 162 ITER=3 l: 163 ITER=2 l: 164 ITER=2 l: 165
ITER=2 l: 166 ITER=3 l: 167 ITER=2 l: 168 ITER=1 l: 169 ITER=1
l: 170 ITER=0 l: 171 ITER=3 l: 172 ITER=1 l: 173 ITER=2 l: 174
ITER=1 l: 175 ITER=1 l: 176 ITER=0 l: 177 ITER=3 l: 178 ITER=2
l: 179 ITER=4 l: 180 ITER=2 l: 181 ITER=2 l: 182 ITER=0 l: 183
ITER=3 l: 184 ITER=3 l: 185 ITER=1 l: 186 ITER=1 l: 187 ITER=1
l: 188 ITER=2 l: 189 ITER=1 l: 190 ITER=2 l: 191 ITER=1 l: 192
ITER=2 l: 193 ITER=1 l: 194 ITER=2 l: 195 ITER=1 l: 196 ITER=3
l: 197 ITER=1 l: 198 ITER=2 l: 199 ITER=1 l: 200 ITER=3 l: 201
ITER=1 l: 202 ITER=0 l: 203 ITER=0 l: 204 ITER=3 l: 205 ITER=2
l: 206 ITER=1 l: 207 ITER=0 l: 208 ITER=1 l: 209 ITER=1 l: 210
ITER=0 l: 211 ITER=2 l: 212 ITER=2 l: 213 ITER=2 l: 214 ITER=1
l: 215 ITER=2 l: 216 ITER=2 l: 217 ITER=2 l: 218 ITER=2 l: 219
ITER=0 l: 220 ITER=3 l: 221 ITER=1 l: 222 ITER=0
===
Subject:
Re: Linking Clapack library in Visual.netHave you tried
downloading CLAPACK3-Windows.zip from www.netlib.org? Itmight
have information on the proper link settings (for VS6
anyway,hopefully it will translate successfully to VS.NET).> I
am visual.net and I successfully compiled the lapack library.>
However when I introduce it to my code I have some library
conflict> such as:> Linking...> MSVCRTD.lib(MSVCR71D.dll) :
error LNK2005: _printf already defined in>
LIBCD.lib(printf.obj)> MSVCRTD.lib(MSVCR71D.dll) : error
LNK2005: _free already defined in> LIBCD.lib(dbgheap.obj)>
MSVCRTD.lib(MSVCR71D.dll) : error LNK2005: _fprintf already
defined in> LIBCD.lib(fprintf.obj)> MSVCRTD.lib(MSVCR71D.dll)
: error LNK2005: _malloc already defined in>
LIBCD.lib(dbgheap.obj)> MSVCRTD.lib(MSVCR71D.dll) : error
LNK2005: _exit already defined in> LIBCD.lib(crt0dat.obj)>
MSVCRTD.lib(ti_inst.obj) : error LNK2005: private: __thiscall>
type_info::type_info(class type_info const &)>
(??0type_info@@AAE@ABV0@@Z) already defined in
LIBCD.lib(typinfo.obj)> MSVCRTD.lib(ti_inst.obj) : error
LNK2005: private: class type_info &> __thiscall
type_info::operator=(class type_info const &)>
(??4type_info@@AAEAAV0@ABV0@@Z) already defined in>
LIBCD.lib(typinfo.obj)> LINK : warning LNK4098: defaultlib
'MSVCRTD' conflicts with use of> other libs; use
/NODEFAULTLIB:library> Debug/PP_stat.exe : fatal error
LNK1169: one or more multiply defined> symbols found> Do you
why is that?> What should I do?> Regards> David
simonin
===
Subject: Re: Linking Clapack library in Visual.netI
think he compiled the main project and the library
withdifferent compiler options. MSVCRTD is for dll possibly
dueto MFC and libcd is for static library.> Have you tried
downloading CLAPACK3-Windows.zip from www.netlib.org? It>
might have information on the proper link settings (for VS6
anyway,> hopefully it will translate successfully to
VS.NET).>>I am visual.net and I successfully compiled the
lapack library.However when I introduce it to my code I have
some library conflict> such
as:>>Linking...MSVCRTD.lib(MSVCR71D.dll) : error LNK2005:
_printf already defined in>
LIBCD.lib(printf.obj)MSVCRTD.lib(MSVCR71D.dll) : error
LNK2005: _free already defined in>
LIBCD.lib(dbgheap.obj)MSVCRTD.lib(MSVCR71D.dll) : error
LNK2005: _fprintf already defined in>
LIBCD.lib(fprintf.obj)MSVCRTD.lib(MSVCR71D.dll) : error
LNK2005: _malloc already defined in>
LIBCD.lib(dbgheap.obj)MSVCRTD.lib(MSVCR71D.dll) : error
LNK2005: _exit already defined in>
LIBCD.lib(crt0dat.obj)MSVCRTD.lib(ti_inst.obj) : error
LNK2005: private: __thiscall> type_info::type_info(class
type_info const &)> (??0type_info@@AAE@ABV0@@Z) already
defined in LIBCD.lib(typinfo.obj)MSVCRTD.lib(ti_inst.obj) :
error LNK2005: private: class type_info &> __thiscall
type_info::operator=(class type_info const &)>
(??4type_info@@AAEAAV0@ABV0@@Z) already defined in>
LIBCD.lib(typinfo.obj)LINK : warning LNK4098: defaultlib
'MSVCRTD' conflicts with use of> other libs; use
/NODEFAULTLIB:libraryDebug/PP_stat.exe : fatal error LNK1169:
one or more multiply defined> symbols found>>Do you why is
that?What should I do?>>Regards>>David simonin>
===
Subject:
Using QR in LAPACKHi allI am struggling somewhat to reproduce
the effects of the qr decompositionfunction in MATLAB with the
Intel MKL LAPACK library. I don't think myquestions are
specific to the Intel variant, however.I am using the
following 3x3 test matrix: 1.0000 0.5000 - 0.3000i 0.3333 -
0.2000i -0.5000 - 0.1000i 1.1000 0.6667 - 0.0500i 0.3333
0.6667 - 0.1000i 1.0000MATLAB gives the QR decomposition as:Q
= -0.8540 -0.2742 + 0.0805i -0.4330 + 0.0386i 0.4270 + 0.0854i
-0.8465 - 0.0059i -0.3054 + 0.0200i -0.2846 -0.4490 + 0.0034i
0.8469 + 0.0059iR = -1.1709 -0.1471 + 0.1907i -0.2889 +
0.0925i 0 -1.3922 + 0.0911i -1.1206 + 0.0708i 0 0 0.4903+
0.0698iA quick multiplication of Q and R gives the expec
product.A subset of my code: int info = 0; int lda = N1; int
lwork = M1*N1; complex<FTYPE> * pTau = new
complex<FTYPE>[M1*N1]; complex<FTYPE> * pWork = new
complex<FTYPE>[M1*N1*16]; cgeqrf(&M1, &N1, temp.m_cfArray,
&lda, pTau, pWork, &lwork, &info); cout << lock << info = <<
info << endl << unlock;// cungrq(&M1, &N1, &M1,
temp.m_cfArray, &lda, pTau, pWork, &lwork, &info);// cout <<
lock << lwork = << lwork << endl << unlock;// cout << lock <<
info = << info << endl << unlock;The result (temp.m_cfArray)
immediately following the call to cgeqrf is asfollows (column
major ordering assumed): -1.1709 - 0.0000i -0.1471 + 0.1907i
-0.2889 + 0.0925i -0.2303 - 0.0461i -1.3952 - 0.0000i -1.1229
- 0.0025i 0.1535 - 0.0000i 0.2128 - 0.0026i -0.4952 -
0.0000iThe LAPACK docs indicate that R is returned in the
upper triangle. This iscorrect for the first row only.The docs
also indicate that the cungrq function is needed to generate
the Qmatrix from the elementary reflectors stored in the lower
triangle of theresult. I have been unable to make this work.I
would greatly appreciate it if someone could show me how to
obtain thesame results MATLAB does. This should be fairly
straightforward, I imagine -but so far, no luck.Thanks in
advance for any and all replies!Ryan
===
Subject: Re: Using QR
in LAPACKRyan, I have no problem reproducing the Matlab
results. Suggestions: 1. Make sure you are calling CGEQRF to
get the QR factorization 2. Make sure you are calling CUNGQR
to get back the Q matrix (not CUNGRQas in your original post).
3. Make sure you are not messing up the 2-d array conventions
between Cand Fortran: Fortran is column major order, C is row
major order. In any case, here are my results. For your test
matrix in Matlab I get(note that the number of significant
digits may be different than yours):m = 1.0000 0.5000 -
0.3000i 0.3333 - 0.2000i -0.5000 - 0.1000i 1.1000 0.6667 -
0.0500i 0.3333 0.6667 - 0.1000i 1.0000q = -0.8540 -0.2684 +
0.0983i 0.4341 + 0.0228i 0.4270 + 0.0854i -0.8451 + 0.0494i
0.3052 + 0.0233i -0.2846 -0.4479 + 0.0328i -0.8376 - 0.1252ir
= -1.1709 -0.1471 + 0.1907i -0.2889 + 0.0925i 0 -1.3952
-1.1229 - 0.0025i 0 0 -0.4952 My Fortran 95 test program
gives:m = 1.0000 + 0.0000i 0.5000 + -0.3000i 0.3333 + -0.2000i
-0.5000 + -0.1000i 1.1000 + 0.0000i 0.6667 + -0.0500i 0.3333 +
0.0000i 0.6667 + -0.1000i 1.0000 + 0.0000iq= -0.8540 + 0.0000i
-0.2684 + 0.0983i 0.4341 + 0.0228i 0.4270 + 0.0854i -0.8451 +
0.0494i 0.3052 + 0.0233i -0.2846 + 0.0000i -0.4479 + 0.0328i
-0.8376 + -0.1252ir= -1.1709 + 0.0000i -0.1471 + 0.1907i
-0.2889 + 0.0925i 0.0000 + 0.0000i -1.3952 + 0.0000i -1.1229 +
-0.0025i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.4952 +
0.0000iwhich is identical to the Matlab output to 4 digits
(all that I prin fromFortran).Here's my Fortran code:program
testqr implicit none complex(4) :: m(3, 3) complex(4) ::
tau(3) complex(4) :: work(3) complex(4) :: r(3, 3) = (0.0,
0.0) integer(4) :: info, ii, jj! test matrix m(1, 1) =
(1.0000, 0.0000) m(1, 2) = (0.5000, -0.3000) m(1, 3) =
(0.3333, -0.2000) m(2, 1) = (-0.5000, -0.1000) m(2, 2) =
(1.1000, 0.0000) m(2, 3) = (0.6667, -0.0500) m(3, 1) =
(0.3333, 0.0000) m(3, 2) = (0.6667, -0.1000) m(3, 3) =
(1.0000, 0.0000) call printcomplex(3, 3, m) print *! lapack
3.0 qr decomposition call CGEQRF( 3, 3, m, 3, TAU, WORK, 3,
INFO )! extract upper triangular part of m to get r do jj = 1,
3 do ii = 1, jj r(ii, jj) = m(ii, jj) end do end do call
printcomplex(3, 3, r) print *! get q matrix call CUNGQR(3, 3,
3, m, 3, TAU, WORK, 3, INFO ) call printcomplex(3, 3, m) print
* contains subroutine printcomplex(m, n, a) integer(4),
intent(in) :: m integer(4), intent(in) :: n complex(4),
intent(in) :: a(m, n) integer(4) :: ii, jj do ii = 1, 3 print
'(3(f8.4,a,f8.4,a))', (real(a(ii,jj)), + , & aimag(a(ii,jj)),
i , jj=1,3) end do end subroutine printcomplexend program
testqrGood luck,John> Hi all> I am struggling somewhat to
reproduce the effects of the qr decomposition> function in
MATLAB with the Intel MKL LAPACK library. I don't think my>
questions are specific to the Intel variant, however.> I am
using the following 3x3 test matrix:> 1.0000 0.5000 - 0.3000i
0.3333 - 0.2000i> -0.5000 - 0.1000i 1.1000 0.6667 - 0.0500i>
0.3333 0.6667 - 0.1000i 1.0000> MATLAB gives the QR
decomposition as:> Q => -0.8540 -0.2742 + 0.0805i -0.4330 +
0.0386i> 0.4270 + 0.0854i -0.8465 - 0.0059i -0.3054 + 0.0200i>
-0.2846 -0.4490 + 0.0034i 0.8469 + 0.0059i> R => -1.1709
-0.1471 + 0.1907i -0.2889 + 0.0925i> 0 -1.3922 + 0.0911i
-1.1206 +0.0708i> 0 00.4903> + 0.0698i> A quick multiplication
of Q and R gives the expec product.> A subset of my code:> int
info = 0;> int lda = N1;> int lwork = M1*N1;> complex<FTYPE> *
pTau = new complex<FTYPE>[M1*N1];> complex<FTYPE> * pWork = new
complex<FTYPE>[M1*N1*16];> cgeqrf(&M1, &N1, temp.m_cfArray,
&lda, pTau, pWork, &lwork, &info);> cout << lock << info = <<
info << endl << unlock;> // cungrq(&M1, &N1, &M1,
temp.m_cfArray, &lda, pTau, pWork, &lwork,&info);> // cout <<
lock << lwork = << lwork << endl << unlock;> // cout << lock
<< info = << info << endl << unlock;> The result
(temp.m_cfArray) immediately following the call to cgeqrf is
as> follows (column major ordering assumed):> -1.1709 -
0.0000i -0.1471 + 0.1907i -0.2889 + 0.0925i> -0.2303 - 0.0461i
-1.3952 - 0.0000i -1.1229 - 0.0025i> 0.1535 - 0.0000i 0.2128 -
0.0026i -0.4952 - 0.0000i> The LAPACK docs indicate that R is
returned in the upper triangle. This is> correct for the first
row only.> The docs also indicate that the cungrq function is
needed to generate theQ> matrix from the elementary reflectors
stored in the lower triangle of the> result. I have been unable
to make this work.> I would greatly appreciate it if someone
could show me how to obtain the> same results MATLAB does.
This should be fairly straightforward, Iimagine -> but so far,
no luck.> Thanks in advance for any and all replies!>
Ryan
===
Subject: Re: Compiling Lapack Blas with Intel Fortran
(IFC)> MKL is free? I thought it was $199.00.> Oh sorry, I was
referring to the academic license which I believe is still
free.I searched on the intel website, but I could not find a
free academiclicense version for the MKL. Where would I find
this, or has it beendiscontinued?Ushnish
===
Subject:
Eigenvalues of reduced matrixfind any solution for it.The
basic question i have is, how do eigenvalues of a reduced
matrixrelate to the eigenvalues of the original matrix. Here
is an examplefor this in Matlab code:M1 = rand(4);EigenM1 =
eig(M1);What surpises me now is, when i leave out one row and
one column thatthe new eigenvalues are completely different to
the ones before. ThatisEigenM2 = eig(M1([1 2 4],[1 2 4]));Is
there any explanation how values of EigenM1 relate to
EigenM2?Thank you and have a nice day.Vedran
===
Subject: Re:
Eigenvalues of reduced matrixhello again and thank you for
your replies.arnold's reference to horn and johnson is
helpfull, but i still look for thesolution of my problem. what
i have found in the book is the followingtheorem.*********Let A
be a hermitian matrix, let r be an integer with 1<=r<=n, and
let A_rdenote any r-by-r principal submatrix of A (obtained by
deleting n-r rowsand the corresponding columns from A). For
each integer k such that 1<=k<=rwe haveLambda_k(A) <=
Lambda_k(A_r) <= Lambda_k+n-r(A)*********I still see no way
how to calculate eigenvalues of the principal submatrix,given
the eigenvalues of the full (in this case symmetrical) matrix.
I needto perform this calculation in order to verify an
implementation of abilinear approximation for a non-linear
system.This still looks like a mathematical challenge to me.
Do you have anyopinions on this?Vedran
===
Subject: Re:
Eigenvalues of reduced matrix> Let A be a hermitian matrix,
let r be an integer with 1<=r<=n, and let A_r> denote any
r-by-r principal submatrix of A (obtained by deleting n-r
rows> and the corresponding columns from A). For each integer
k such that 1<=k<=r> we have> Lambda_k(A) <= Lambda_k(A_r) <=
Lambda_k+n-r(A)Yes, that's all one can say in general. There
are no other relations,not even in the 2 by 2 case, where it
is easy to convince yourself!> I still see no way how to
calculate eigenvalues of the principal submatrix,> given the
eigenvalues of the full (in this case symmetrical) matrix. I
need> to perform this calculation in order to verify an
implementation of a> bilinear approximation for a non-linear
system.Simply calculate the smaller problem independently of
the larger one.This is almost the best you can do. With more
understanding you could alsosolve the reduced problem as a
rank 1 update problem, which would save thecomputer a little
effort at the expense of much more programming,wasting much
more of your time...(See Section 8.6.3 in the book by Golub
and van Loan)Arnold Neumaier
===
Subject: Re: Eigenvalues of
reduced matrix> find any solution for it.> The basic question
i have is, how do eigenvalues of a reduced matrix> relate to
the eigenvalues of the original matrix. Here is an example>
for this in Matlab code:> M1 = rand(4);> EigenM1 = eig(M1);>
What surpises me now is, when i leave out one row and one
column that> the new eigenvalues are completely different to
the ones before. That> is> EigenM2 = eig(M1([1 2 4],[1 2
4]));> Is there any explanation how values of EigenM1 relate
to EigenM2?> Thank you and have a nice day.> VedranThe
determinant of a general square matrix (an operator in some
finitebasis) is the volume of a parallelepiped whose edges are
the columns of thematrix. When you delete a row and its
corresponding column the determinantof the reduced matrix will
differ dimensionally if not also numericallyfrom the original.
If the matrix is self adjoint then the interleavingtheorem of
Hylleraas, Undheim, and MacDonald (ca 1930's in atomic
structurecontext where the name of the game is to minimize the
volume by guessing abasis that is as complete as one can afford
computationally) says that thevolume will shrink if you expand
the matrix by the addition of a row andcolumn.-- HTH,Gerry
T.
===
Subject: Re: Eigenvalues of reduced matrix> The basic
question i have is, how do eigenvalues of a reduced matrix>
relate to the eigenvalues of the original matrix. Here is an
example> for this in Matlab code:> M1 = rand(4);> EigenM1 =
eig(M1);> What surpises me now is, when i leave out one row
and one column that> the new eigenvalues are completely
different to the ones before.A physical analogy is handy if
your matrix is real symmetric or can bewritten as the product
of two real symmetric matrices (the generalizedeigenvalue
problem). The matrix can then be thought of as a
discreteelement vibration model of some mechanical system
(e.g., springs andmasses), with the eigenvalues proportional
to the square of the naturalfrequencies. (Positive
definiteness also helps in this analogy.) DeletingRow i and
Column i of the matrix corresponds physically to applying
aconstraint to DOF i, thus changing the system and its
eigenvalues.
===
Subject: Re: Eigenvalues of reduced matrix>
find any solution for it.> The basic question i have is, how
do eigenvalues of a reduced matrix> relate to the eigenvalues
of the original matrix. Here is an example> for this in Matlab
code:> M1 = rand(4);> EigenM1 = eig(M1);> What surpises me now
is, when i leave out one row and one column that> the new
eigenvalues are completely different to the ones before. That>
is> EigenM2 = eig(M1([1 2 4],[1 2 4]));> Is there any
explanation how values of EigenM1 relate to EigenM2?For
symmetric matrices, there are well-known interlacing
results(which you can find, e.g. in Parlett, or in Horn and
Johnson).In the nonsymmetric case, anything can happen.Arnold
Neumaier
===
Subject: Re: AT LAST - TRADING THE MARKETS BY
UNIFYING CYCLES AND PROBABILITY!!!!send in fo
===
Subject:
Indexes of uique product termshere is a problem of how to
determine coefficients containing unique valuesin a vector
resulting by taking a Kronecker product of a vector.
Thefollowing examples in matlab notation explains this.>> syms
x1 x2>> x=[x1;x2];>> y=kron(x,x)y =[ x1^2][ x1*x2][ x1*x2][
x2^2]obviously only y(1), y(2) and y(4) represent unique terms
of this product.( ix = [1,2,4] )Also>> syms x1 x2 x3>>
x=[x1;x2;x3];>> y=kron(x,x)y =[ x1^2][ x1*x2][ x1*x3][ x1*x2][
x2^2][ x2*x3][ x1*x3][ x2*x3][ x3^2]The relevant idices here
are ix = [1,2,3,5,6,9] and so on.One way to find these is to
think of y as a 3 x 3 (length(x) x length(x))matrix and to
take a lower triangular matrix to determine indices of
uniqueproduct terms. This even works for multiple kronecker
product such as,kron(x,kron(x,x)) if you think of the
triangular matrix in higherdimensions.My question now is, does
anybody know any other, possibly more efficient wayto obtain
these indices. Your comments are welcome.Have a nice dayVedran
Dizdarevic
===
Subject: Announce: CamlFloatWe are glad to
announce the first release of CamlFloat, a highlystylized
OCaml interface to Lapack and Blas. It can be used
eitherinteractively (like Matlab) or as a library in your own
programs. Itcomes with a tutorial and plenty of
documentation.OCaml (http://www.ocaml.org/) is a functional,
imperative,object-orien, strict, static type-inferring
language in the MLfamily. It compiles both to byte-code and
native code on almostall modern processors and operating
systems.CamlFloat tries to create a Matlab-like environment in
OCaml withoutany attendant loss of speed or memory efficiency.
And of course withfull access to all the power of OCaml.The
library has been used extensively to code, test and time our
ownresearch algorithms. So it is battle-hardened to some
extent. Userfeedback is welcome.The package can be obtained
from http://math.ucsb.edu/~lyons/camlFloat/Enjoy,Shivkumar
ChandrasekaranBill Lyons