//
//<![CDATA[
//
//  Algoritmo di tridiagonalizzazione di una matrice arbitraria.
//  Scopritore ( o inventore ) Giampaolo Bottoni.
//  Implementazione del 31 agosto 2008.
//
//  Se la tridiagonalizzazione non riesce il risultato è comunque
//  una matrice in forma di Hessenberg.
//  Dunque IN OGNI CASO la matrice trasformata e' adatta per venire 
//  trattata con l'algoritmo QR di Francis per la ricerca di autovalori 
//  e autovettori.
//
//
var rinuncio_se=1.e-9;
var mifermoqui=1000000;
//
function tridiago(a){
   //
   // a deve essere un vettore di vettori che in uscita conterra'
   // la tridiagonale non simmetrica.
   // 
   var j,k,p,t,s="",n=a.length;
   var x=new Array(n);
   var y=new Array(n);
   var z=new Array(n);
   //
   for(p=1;n-1>p;p++){
       //
       // Questo serve per vedere cosa è successo dopo un certo numero p di passi...
       //
       if(p>mifermoqui) break;
       //
       // Opera sulla colonna p-1 _esima.
       //
       for(j=p;n>j;j++){
           x[j-p]=a[j][p-1];
           }
       for(k=0;n>k;k++){
           for(j=p;n>j;j++){
               y[j-p]=a[j][k];
               }
           householder(z,x,y,n-p); // function definita poi...
           for(j=p;n>j;j++){
               a[j][k]=z[j-p];
               if(1.0e-14>Math.abs(z[j-p])) a[j][k]=0;
               }
           }
       for(k=0;n>k;k++){
           for(j=p;n>j;j++){
               y[j-p]=a[k][j];
               }
           householder(z,x,y,n-p);
           for(j=p;n>j;j++){
               a[k][j]=z[j-p];
               if(1.0e-14>Math.abs(z[j-p])) a[k][j]=0;
               }
           }
       //
       // Vede se la riga è gia nulla. Prosegue ad azzerarla
       // solo se c'è veramente bisogno.
       //
       t=0.0;
       for(j=p+1;n>j;j++)t+=Math.abs(a[p-1][j]);
       if(1.e-14*Math.abs(a[p][p-1])>t) continue; 
       //
       // Opera sulla riga - parte aggiuntiva rispetto alla 
       // normale trasformazione in forma di Hessenberg.
       //
       if(n>p+2){
           for(j=p+1;n>j;j++){
               x[j-p-1]=a[p-1][j];
               }
           for(k=0;n>k;k++){
               for(j=p+1;n>j;j++){
                   y[j-p-1]=a[j][k];
                   }
               householder(z,x,y,n-p-1);
               for(j=p+1;n>j;j++){
                   a[j][k]=z[j-p-1];
                   if(1.0e-14>Math.abs(z[j-p-1])) a[j][k]=0;
                   }
               }
           for(k=0;n>k;k++){
               for(j=p+1;n>j;j++){
                   y[j-p-1]=a[k][j];
                   }
               householder(z,x,y,n-p-1);
               for(j=p+1;n>j;j++){
                   a[k][j]=z[j-p-1];
                   if(1.0e-14>Math.abs(z[j-p-1])) a[k][j]=0;
                   }
               }
           }
       //
       // Dovrei usare a[p-1][p] a denominatore. Per evitare
       // valori troppo grandi potrei rinunciare a tridiagonalizzare
       // quando a[p-1][p] è troppo piccolo... ma invece...
       //
       if(1.e-30 > Math.abs(a[p-1][p])) a[p-1][p]=1.e-30;
       //
       // ATTENZIONE: in base alla sperimentazione numerica
       // sembra corretto sostenere che gli autovalori della
       // tridiagonale sono molto vicini a quelli della matrice di
       // partenza ANCHE SE risulta che t è praticamente infinito. 
       // In base a queste osservazioni ho aggiunto il solo "vincolo" 
       // che l'elemento a[p-1][p] non possa essere esattamente zero
       // al solo scopo di limitare empiricamente i valori assumibili
       // da t.
       //
       // In sostanza mi sembra di dover desumere che, nel caso in cui
       // l'elemento a[p-1][p] sia zero, la tridiagonale è PARTIZIONABILE
       // il che dovrebbe rendere questa situazione molto auspicabile
       // e assolutamente non preoccupante.
       //
       if(Math.abs(a[p-1][p])>rinuncio_se){
           //
           // Qui eseguo la trasformazione "delicata".
           // Se salto questa parte viene fuori una normale
           // trasformazione in forma di Hessenberg ma
           // appesantita dall'aver puntato all'azzeramento
           // delle righe oltre che delle colonne.
           //
           // I valori cruciali sono:
           // alert(a[p-1][p-1]+", "+a[p-1][p]+", "+a[p-1][p+1]+"\n"+a[p][p-1]);
           //
           t=a[p-1][p+1]/a[p-1][p];
           for(k=0;n>k;k++){
               a[p][k]+=t*a[p+1][k];
               }
           for(k=0;n>k;k++){
               a[k][p+1]-=t*a[k][p];
               }
           // 
           //  Questo azzeramento sarebbe inutile ma per ragioni numeriche
           //  diventa necessario per evitare che valori trascurabili possano
           //  però creare delle difficoltà nei passi successivi...
           //
           for(k=p+1;n>k;k++){
               a[p-1][k]=0;
               a[k][p-1]=0;
               }
           // 
           // Memorizzo la storia delle trasformazioni...
           //
           s+="Al passo "+p+") t = "+t+
              String.fromCharCode(60)+"br/"+String.fromCharCode(62);
           }
       else {
           s+="Al passo "+p+") rinuncia a tridiagonalizzare"+
              String.fromCharCode(60)+"br/"+String.fromCharCode(62);
           }
       }  // fine del passo p_esimo.
   return s;
   }
function hopre(x){
   //
   // Prepara il vettore x[] in modo che sia usabile
   // come trasformatore della matrice di Householder che
   // moltiplicata per x[] fornisce come risultato la 
   // norma di x[] nella componente zeresima e rende nulle
   // tutte le componenti esclusa quella iniziale. 
   // In x[0] (se non c'e') mette un piccolo vettore
   // contenente i dati per fare la trasformazione ossia
   // in x[0][0] la vera componente iniziale del vettore x[],
   // in x[0][1] il fattore moltiplicativo, in x[0][2]
   // la componente iniziale del vettore della diade e in
   // x[0][3] il segno di x[0][0]. Pertanto se x[0][3] è
   // nullo vuol dire che la preparazione va rifatta. 
   // 
   var m,r,j,xx,sg;
   m=x.length;
   if(hopre.arguments.length>1) m=hopre.arguments[1];
   if(isNaN(x[0])) { xx=x[0][0];}
   else { xx=x[0]; x[0]=[xx,0,0,0];}
   r=xx*xx;
   for(j=1;m>j;j++){
       r+=x[j]*x[j]; }
   r=Math.sqrt(r);
   if(0>=r) { xx=1.e-50; r=xx; }
   sg=1;
   if(0>xx) sg=-1;           
   x[0][1]=sg/(r*(r+Math.abs(xx)));
   x[0][2]=xx+sg*r;
   x[0][3]=sg;
   return xx;
   }
//
function householder(z,x,y){ // Risultato nel vettore z[]
   var j,k,n,sc;
   n=x.length;
   if(householder.arguments.length>3) n=householder.arguments[3];
   // Quando si vuole usare una nuova trasformazione
   // la componente x[0] deve essere un numero o, se 
   // fosse un vettore, bisogna che x[0][3] sia nullo.
   if(!isNaN(x[0])) hopre(x,n);
   else if( 0.5 > Math.abs(x[0][3]) ) hopre(x,n);
   sc=x[0][2]*y[0];
   for(j=1;n>j;j++) sc+=x[j]*y[j];
   sc*=x[0][1];
   if(0>x[0][0]){
       z[0]=x[0][2]*sc +y[0];
       for(k=1;n>k;k++){
           z[k]=x[k]*sc+y[k];        
           }
       }
   else {
       z[0]=x[0][2]*sc-y[0];
       for(k=1;n>k;k++){
           z[k]=x[k]*sc-y[k];        
           }
       }
   return n;
   }
//
// ]]> Fine dell'algoritmo di tridiagonalizzazione
//