Index: Singular/LIB/jordan.lib
===================================================================
--- Singular/LIB/jordan.lib	(revision d694de170d576d5f00d6f1047b53246e20433998)
+++ Singular/LIB/jordan.lib	(revision 092430050dbf6ad127e0cf157b408a5b79bf6403)
@@ -1,5 +1,5 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-version="$Id: jordan.lib,v 1.11 1999-06-16 17:24:35 mschulze Exp $";
+version="$Id: jordan.lib,v 1.12 1999-06-21 08:45:41 mschulze Exp $";
 info="
 LIBRARY: jordan.lib  PROCEDURES TO COMPUTE THE JORDAN NORMAL FORM
@@ -17,14 +17,14 @@
 static proc countblocks(matrix M)
 {
-  int n;
+  int b,r,r0;
 
   int i=1;
-  int r,r0;
   while(i<=nrows(M))
   {
-    n++;
+    b++;
     r=nrows(M[i]);
     r0=r;
 
+    dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
     while(i<r0&&i<nrows(M))
     {
@@ -39,9 +39,10 @@
       }
     }
+    dbprint(printlevel-voice+2,"//...block "+string(b)+" found");
 
     i++;
   }
 
-  return(n);
+  return(b);
 }
 ///////////////////////////////////////////////////////////////////////////////
@@ -72,5 +73,4 @@
 "
 {
-  // check parameter M
   int n=nrows(M);
   if(n==0)
@@ -85,28 +85,54 @@
   }
 
-  // replace M by its constant part
   M=jet(M,0);
 
-  // try to increase number of blocks by transposing M
-  if(countblocks(M)<countblocks(transpose(M)))
-  {
+  dbprint(printlevel-voice+2,"//counting blocks of matrix...");
+  int i=countblocks(M);
+  dbprint(printlevel-voice+2,"//...blocks of matrix counted");
+  if(i==1)
+  {
+    dbprint(printlevel-voice+2,"//matrix has 1 block");
+  }
+  else
+  {
+    dbprint(printlevel-voice+2,"//matrix has "+string(i)+" blocks");
+  }
+
+  dbprint(printlevel-voice+2,"//counting blocks of transposed matrix...");
+  int j=countblocks(transpose(M));
+  dbprint(printlevel-voice+2,"//...blocks of transposed matrix counted");
+  if(j==1)
+  {
+    dbprint(printlevel-voice+2,"//transposed matrix has 1 block");
+  }
+  else
+  {
+    dbprint(printlevel-voice+2,"//transposed matrix has "+string(j)+" blocks");
+  }
+
+  if(i<j)
+  {
+    dbprint(printlevel-voice+2,"//transposing matrix...");
     M=transpose(M);
-  }
-
-  // compute eigenvalues eM and multiplicities mM of M
-  list fac;
+    dbprint(printlevel-voice+2,"//...matrix transposed");
+  }
+
+  list fd;
   matrix M0;
+  poly cp;
   ideal eM,eM0;
   intvec mM,mM0;
   intvec u;
-  int i=1;
-  int j,r,r0;
+  int b,r,r0;
+
+  i=1;
   while(i<=nrows(M))
   {
+    b++;
     u=i;
     r=nrows(M[i]);
     r0=r;
 
-    // find next block
+    dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
     while(i<r0&&i<nrows(M))
     {
@@ -122,8 +148,11 @@
       }
     }
-
-    // if 1x1-block
+    dbprint(printlevel-voice+2,"//...block "+string(b)+" found");
+
     if(size(u)==1)
     {
+      dbprint(printlevel-voice+2,"//1x1-block:");
+      dbprint(printlevel-voice+2,M[u[1]][u[1]]);
+
       if(mM[1]==0)
       {
@@ -139,10 +168,19 @@
     else
     {
-      // factorize characteristic polynomial of block M0
+      dbprint(printlevel-voice+2,
+        "//"+string(size(u))+"x"+string(size(u))+"-block:");
       M0=getblock(M,u);
-      fac=factorize(det(M0-var(1)*freemodule(size(u))),2);
-      eM0,mM0=fac[1..2];
-
-      // compute eigenvalues eM0 of M0
+      dbprint(printlevel-voice+2,M0);
+
+      dbprint(printlevel-voice+2,"//characteristic polynomial:");
+      cp=det(module(M0-var(1)*freemodule(size(u))));
+      dbprint(printlevel-voice+2,cp);
+
+      dbprint(printlevel-voice+2,"//factorizing characteristic polynomial...");
+      fd=factorize(cp,2);
+      dbprint(printlevel-voice+2,"//...characteristic polynomial factorized");
+
+      dbprint(printlevel-voice+2,"//computing eigenvalues...");
+      eM0,mM0=fd[1..2];
       if(1<var(1))
       {
@@ -169,4 +207,5 @@
         }
       }
+      dbprint(printlevel-voice+2,"//...eigenvalues computed");
 
       if(mM[1]==0)
@@ -185,5 +224,5 @@
   }
 
-  // sort eigenvalues
+  dbprint(printlevel-voice+2,"//sorting eigenvalues...");
   poly e;
   int m;
@@ -203,6 +242,7 @@
     }
   }
-
-  // remove multiple eigenvalues
+  dbprint(printlevel-voice+2,"//...eigenvalues sorted");
+
+  dbprint(printlevel-voice+2,"//removing multiple eigenvalues...");
   i=1;
   j=2;
@@ -217,4 +257,5 @@
       i++;
       eM[i]=eM[j];
+      mM[i]=mM[j];
     }
     j++;
@@ -222,6 +263,11 @@
   eM=eM[1..i];
   mM=mM[1..i];
-
-  // get optional parameter opt
+  dbprint(printlevel-voice+2,"//...multiple eigenvalues removed");
+
+  dbprint(printlevel-voice+2,"//eigenvalues:");
+  dbprint(printlevel-voice+2,eM);
+  dbprint(printlevel-voice+2,"//multiplicities:");
+  dbprint(printlevel-voice+2,mM);
+
   int opt=0;
   if(size(#)>0)
@@ -236,8 +282,6 @@
     return(list(eM));
   }
-
-  // define variables
   int k,l;
-  matrix E=freemodule(n);
+  matrix I=freemodule(n);
   matrix Mi,Ni;
   module sNi;
@@ -256,7 +300,9 @@
   for(i=ncols(eM);i>=1;i--)
   {
-    Mi=M-eM[i]*E;
-
-    // compute kernels K of powers of Mi
+    Mi=M-eM[i]*I;
+
+    dbprint(printlevel-voice+2,
+      "//computing kernels of powers of matrix minus eigenvalue "
+      +string(eM[i]));
     K=list(module());
     for(Ni,sNi=Mi,0;size(sNi)<mM[i];Ni=Ni*Mi)
@@ -265,8 +311,11 @@
       K=K+list(sNi);
     }
+    dbprint(printlevel-voice+2,"//...kernels computed");
 
     if(opt<=1)
     {
-      // compute Jordan block size vector bMi corresponding to eigenvalue eM[i]
+      dbprint(printlevel-voice+2,
+        "//computing Jordan block sizes for eigenvalue "
+        +string(eM[i])+"...");
       bMi=0;
       bMi[size(K[2])]=0;
@@ -279,9 +328,12 @@
       }
       bM=list(bMi)+bM;
+      dbprint(printlevel-voice+2,"//...Jordan block sizes computed");
     }
 
     if(opt>=1)
     {
-      // compute Jordan basis vectors K corresponding to eigenvalue eM[i]
+      dbprint(printlevel-voice+2,
+        "//computing Jordan basis vectors for eigenvalue "
+        +string(eM[i])+"...");
       if(size(K)>1)
       {
@@ -306,4 +358,5 @@
         }
       }
+      dbprint(printlevel-voice+2,"//...Jordan basis vectors computed");
     }
   }
@@ -336,5 +389,4 @@
 "
 {
-  // check parameters
   if(size(jd)<2)
   {
@@ -365,5 +417,4 @@
   }
 
-  // get size of Jordan matrix
   int i,j,k,n;
   for(i=s;i>=1;i--)
@@ -388,5 +439,4 @@
   }
 
-  // create Jordan matrix
   int l;
   matrix J[n][n];
