MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

Please help in creating/installing my package

  • To: mathgroup at smc.vnet.net
  • Subject: [mg79304] Please help in creating/installing my package
  • From: jeremito <jeremit0 at gmail.com>
  • Date: Tue, 24 Jul 2007 06:02:34 -0400 (EDT)

I have searched through this list and on the web, but haven't yet been
able to find an explanation of how to write/install my own packages in
Mathematica.

I have a few functions that I want to be able to access from any
Mathematica document.  I created my package "Arnoldi.m" (copied below)
and put it into my ~/Library/Mathematica/Applications directory.  When
I do

<<Arnoldi.m

There are no errors, but I can't use any of the functions I wrote.
Can someone please help me learn how to write my own packages.
Thanks,
Jeremy

BeginPackage["Arnoldi`"]

arnoldi::usage = "{Q, H} = Arnoldi[A,q,iters]\nwhere A is a matrix, q
\
is the starting vector with same size as A, and iters is the number \
of Arnoldi Iterations.  Arnoldi returns: orthonormal basis vectors, \
columsn of Q; upper Hessenberg matrix, H."

ArnoldiVectors::usage =
 "Vectors = ArnoldiVectors[Q, vecs]\nwhere Q comes an Arnoldi process
\
and vecs are the eigenvectors of H from the same Arnoldi process."

Begin["Arnoldi`"]

arnoldi[A_, v_, iters_] := (
  (* *)
  invariant = False;
  m = Length[A];
  If[iters > m,
   Print["Iterations cannot be greater than size of matrix."];
   iterations = m;
   , iterations = iters];
  (*Print["Iterations = ",iterations]*);

  H = ConstantArray[0, {iterations + 1, iterations}];
  Subscript[q, 1] = v/Norm[v];

  For[k = 1, k <= iterations, k++,
   (*Print["k = ", k];*)
   Subscript[q, k + 1] = A.Subscript[q, k];
   (*Orthogonalize*)
   For[j = 1, j <= k, j++,
    (*Print["j = ",j];*)

    H[[j, k]] = Subscript[q, k + 1].Subscript[q, j];
    Subscript[q, k + 1] =
     Subscript[q, k + 1] - H[[j, k]] Subscript[q, j];
    ];
   (*Normalize*)
   H[[k + 1, k]] = Norm[Subscript[q, k + 1]];

   If[H[[k + 1, k]] < 10^-10,
    Print["I am invariant!"];
    invariant = True;
    Break[],
    Subscript[q, k + 1] = Subscript[q, k + 1]/H[[k + 1, k]];
    ];
   ];
  If[invariant, iterations = k;, iterations = k - 1];
  H = H[[1 ;; iterations + 1, 1 ;; iterations]];
  Q = Table[Subscript[q, i], {i, iterations + 1}];
  Q = Transpose[Q];
  (*vals = Eigenvalues[H[[1;;k-1,1;;k-1]]];*)
  {Q, H}
  )

ArnoldiVectors[Q_, vecs_] := (
  For[j = 1, j <= Length[vecs], j++,
   Subscript[v, j] = ConstantArray[0, Length[vecs[[1]]]];
   For[i = 1, i <= Length[vecs[[j]]], i++,
    Subscript[v, j] = Subscript[v, j] + Q[[All, i]]*vecs[[j, i]];
    ];
   ];
  V1 = Table[Subscript[v, i], {i, Length[vecs]}];
  {V1}
  )

End[]

EndPackage[]



  • Prev by Date: Re: Re: Embedded Style Sheets
  • Next by Date: Re: Mathematica to .NET compiler
  • Previous by thread: Re: a challenging definite integral!
  • Next by thread: RE: Please help in creating/installing my package