Artifact de695296f26c77914fdec4f3d946f905dc356793b4c9475c20b83f2c9f1e7b6c:

• Executable file r37/packages/sparse/spmateig.red — part of check-in [f2fda60abd] at 2011-09-02 18:13:33 on branch master — Some historical releases purely for archival purposes

  git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 7065) [annotate] [blame] [check-ins using] [more...]

• Executable file r38/packages/sparse/spmateig.red — part of check-in [f2fda60abd] at 2011-09-02 18:13:33 on branch master — Some historical releases purely for archival purposes

  git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 7065) [annotate] [blame] [check-ins using]

%***********************************************************************
%=======================================================================
%
% Code for the extension of the Matrix Package to include Sparse 
% Matrices.
%
% The following code is for the functions to calculate eigenvalues and
% the rank of sparse matrices.
%
% Author: Stephen Scowcroft.                           Date: June 1995
%
%=======================================================================
%***********************************************************************

module spmateig;

% Based on the current eigenvalue code.

symbolic procedure spmateigen(u,eival);
%    U is a matrix form, eival an indeterminate naming the eigenvalues.
%    Result is a list of lists:
%      {{eival-eq1,multiplicity1,eigenvector1},....},
%    where eival-eq is a polynomial and eigenvector is a matrix.
%    How much should we attempt to solve the eigenvalue eq.? sqfr?
%    Sqfr is necessary if we want to have the full eigenspace. If there
%    are multiple roots another pass through eigenvector calculation
%    is needed(done).
%    We should actually perform the calculations in the extension
%    field generated by the eigenvalue equation(done inside).
  begin scalar arbvars,exu,sgn,q,r,s,x,y,z,eivec,!*factor,!*sqfree,
               !*exp,val,res,res1,rl;
        integer l,cnt;
     !*exp := t;
     if not(getrtype u eq 'sparse) then typerr(u,"sparse matrix");
     eival := !*a2k eival;
     kord!* := eival . kord!*;
     rl:=sprow_dim(u);
     exu := spmateigen1(spmatsm u,eival);
     q := car exu;
     y := cadr exu;
     z := caddr exu;
     exu := cdddr exu;
     !*sqfree := t;
     for each j in cdr fctrf numr subs2(lc z ./ 1)
          do if null domainp car j and mvar car j eq eival
                then s := (if null red car j
                              then !*k2f mvar car j . (ldeg car j*cdr j)
                            else j) . s;
     for each j in q
       do (if x then rplacd(x,cdr x + cdr j)
             else s := (y . cdr j) . s)
            where x := assoc(y,s) where y := absf reorder car j;
     l := length s;
     r := 'list .
       for each j in s collect
     <<if null((cdr j = 1) and (l = 1)) then
         <<y := 1;
           for each k in exu do
         if x := reduce!-mod!-eig(car j,c!:extmult(k,y))
           then y := x>>;
       arbvars := nil;
       for each k in lpow z do
          if (y=1) or null(k member lpow y) then
         arbvars := (k . makearbcomplex()) . arbvars;
       sgn := (y=1) or evenp length lpow y;
       cnt:=0;
       for each k in lpow z do
        << if x := assoc(k,arbvars) then 
              res:=list(cnt:=cnt+1,(1 . mvar cdr x))
            else 
             << val:= mkgleig(k,y,sgn := not sgn,arbvars);
                if (val=simp 0) then cnt:=cnt+1
                 else res:=list(cnt:=cnt+1,(1 . prepsq!* val)); 
             >>;
          if res=nil then nil 
           else << res1:=append(res1,list(res));
                   res:=nil;
                   eivec:=mkempspmat(rl,list('spm,rl,1));
                     for i:=1:rl do
                      <<letmtr3(list(eivec,i),atsoc(i,res1),eivec,nil);>>;
                >>;
          >>;
        res1:=nil;
       list('list,prepsq!*(car j ./ 1),cdr j,eivec)>>;
     kord!* := cdr kord!*;
     return r
   end;

symbolic procedure spmateigen1(u,eival);
   begin scalar diag,q,x,y,z,w,j,res; integer l,lm,m,cc;
      lm := spcol_dim(u);
      z := 1;
     for rr:=1:sprow_dim(u) do
     << y := 1;
       diag:=nil;
       cc:=findrow(u,rr);
      if not (cc=nil) then
      <<for each xx in cdr cc do
       <<w:=simp cdr xx;
         y := lcm(y,denr w);
       >>;
      >>;
        m := lm;
        l := l + 1;
        x := nil;
       if cc=nil then cc:=list(list(nil),list(nil));
      for each xx in reverse cdr cc do
        <<if xx='(nil) then <<m:=rr+1; j:=simp nil>>
            else << j:=simp cdr xx; m:=car xx>>;
          if not diag then
          <<if m=rr then <<diag:=j; j:=simp nil>>
              else if m<rr then nil
               else <<diag:=simp findelem2(u,rr,rr); m:=rr>>; 
             x:=list m .* multf(addf(numr diag,negf multf(!*k2f eival,
                                denr diag)),quotf(y,denr diag)) .+ x;
             if not (xx=nil) then m:=car xx;
          >>;
          if numr j and not (m=rr) then
            x := list m .* multf(numr j,quotf(y,denr j)) .+ x;
        >>;
        y := z;
        z := c!:extmult(if null red x then <<
               q := (if p then (car p  . (cdr p + 1)) . delete(p,q)
                      else (lc x  . 1) . q) where p = assoc(lc x,q);
                        !*p2f lpow x>> else x,z);
        res:=append(res,list x);
      >>;
       u:=res;
      return q . y . z . u
   end;

flag('(spmateigen),'opfn);

flag('(spmateigen),'noval);

% These functions calculate the rank of a sparse matrix.

% To replace current function.
% Extended for sparse matrices.

symbolic procedure rank!-eval u;
   begin scalar n;
      if cdr u then rerror(matrix,17,"Wrong number of arguments")
       else if getrtype (u := car u) eq 'matrix
         then return rank!-matrix spmatsm u
       else if getrtype (u) eq 'sparse
         then return sprank!-matrix spmatsm u       
       else if null eqcar(u := aeval u,'list) then typerr(u,"matrix")
       else return rank!-matrix
         for each row in cdr u collect
           if not eqcar(row,'list)
               then rerror(matrix,15,"list not in matrix shape")
            else <<row := cdr row;
                   if null n then n := length row
                    else if n neq length row
                     then rerror(matrix,151,"list not in matrix shape");
                   for each j in row collect simp j>>
   end;


put('rank,'psopfn,'rank!-eval);

% The function to calculate the rank.
% Based on the matrix rank function.

symbolic procedure sprank!-matrix u;
 begin scalar x,y,z,w,j,cc; integer m,n;
   z := 1;
   for rr:=1:sprow_dim(u) do
   << y := 1;
      cc:=findrow(u,rr);
      if not (cc=nil) then
      <<for each xx in cdr cc do
        <<w:=simp cdr xx;
           y := lcm(y,denr w);
        >>;
      >>;
      m := 1;
      x := nil;
      if not (cc=nil) then
      <<for each xx in cdr cc do
        <<j:=simp cdr xx;
           m:=car xx;
           if numr j then
            x := list m .* multf(numr j,quotf(y,denr j)) .+ x;
        >>;
      >>;
      if y := c!:extmult(x,z) then <<z := y; n := n + 1>>;
    >>;
     return n
    end;

endmodule;

end;

%***********************************************************************
%=======================================================================
%
% End of Code.
%
%=======================================================================
%***********************************************************************