06/11/2011: Power iteration algorithm to compute dominant eigenvector and eigenvalue.
/************************************************************************/
/* 04/24/2011: Modified. */
/* 04/22/2011: Previouly modified. */
/* 04/22/2011: Created. */
/************************************************************************/
/************************************************************************/
/* Program: power(). */
/* Function: Power iteration to compute eigenvectors and eigenvalues.*/
/************************************************************************/
mata:
struct matrix_struct scalar power(real matrix X,
|string scalar nodots, string scalar left){
struct matrix_struct scalar z /* Return structure. */
real matrix E /* Original or transposed edge (arc) list. */
real vector e0 /* Current iteration eigenvector. */
real vector e1 /* Next iteration eigenvector. */
real scalar i /* Row counter. */
real scalar exit /* Boolean exit indicator. */
real scalar norm /* Vector norm (absolute length). */
real scalar m /* Total node number. */
real scalar k /* Iteration counter. */
string scalar o /* Display scalar. */
/**********************************************************************/
/* Begin program. */
/**********************************************************************/
/* Initialize variables. */
m=max(X[.,(1,2)]); e0=J(m,1,1); e1=J(m,1,0)
/* Transpose elements of edge (arc) list if calculating left eigenvecs*/
if(left!="") E=(X[.,2],X[.,1],X[.,3])
else E=X
/**********************************************************************/
/* Inform user of program starting. */
/**********************************************************************/
o="noisily display "+`""""'; stata(o); displayflush()
o="Power iteration algorithm (iterations)"
o="noisily display "+`"""'+o+`"""'
stata(o); displayflush()
/* Display dots unless specified not to. */
if(nodots==""){
o="----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5"
o="noisily display "+`"""'+o+`"""'
stata(o); displayflush()
}
/* Initialize boolean exit indicator and iteration counters. */
exit=0; k=0
/* Loop to convergence or iteration<=16000 (Mata default). */
while(exit!=1){
++k
/* Loop through row entries of augmented edge list. */
for(i=1;i<=rows(E);i++){
/* E equals augmented edge list. */
/* E[i,1] equals ith entry of A_ij. */
/* E[i,2] equals jth entry of A_ij. */
/* E[i,3] equals ijth entry of A_ij, that is, the weight. */
e1[E[i,1]]=e1[E[i,1]]+E[i,3]*e0[E[i,2]]
}
/***** Calculate next iteration. */
norm=norm(e1); e1=e1/norm
/***** Define exit or continue iteration conditions. */
/* Exit if last and current iteration relatively similar. */
if(mreldif(e0,e1)<(1e-10)|nonmissing(e1)==0){
/* Display status unless specified not to. */
if(nodots=="") display_status(16000,k)
o="noisily display "+`""""'
stata(o); displayflush()
/***** Display failure message if vector values all missing. */
if(nonmissing(e1)==0){
o="Power iteration algorithm failed to converge"
o="noisily display as error "+`"""'+o+`"""'
stata(o); displayflush()
}
/***** Display success message if convergence achieved. */
else{
o="Power iteration algorithm completed"
o="noisily display "+`"""'+o+`"""'
stata(o); displayflush()
}
/****** Return eigenvector and eigenvalue. */
if(nonmissing(e1)==0){
/***** Mata error 3360 - convergence not achieved. */
z.X1=3360; return(z)
}
else{
/***** 1st element is eigenvector, 2nd element is eigenvalue. */
z.X1=e1; z.X2=norm; return(z)
}
exit=1
}
/* If convergence not achieved after 16000 iterations, exit loop. */
else if(k==16000){
/* Display status unless specified not to. */
if(nodots=="") display_status(16000,k)
o="Power iteration algorithm failed to converge"
o="noisily display as error "+`"""'+o+`"""'
stata(o); displayflush()
/***** Mata error 3360 - convergence not achieved. */
z.X1=3360
return(z)
exit=1
}
/* If convergence not achieved, iterate e0 and reset e1. */
else{
e0=e1; e1=J(m,1,0)
/* Display status unless specified not to. */
if(nodots=="") display_status(16000,k)
}
}
}
end