sparse_rand = function( order; density ) {

    # This function returns a random positive-definite matrix with
    # `order' rows and columns and approximately `density*order^2'
    # non-zero elements.

    local( i; v; m );

    v = fill( density*order^2; r*order+0.5 );
    m = 10 * diag( fill( order; r ) );

    for ( i in 1:v.ne ) {
        m[ v[i]; v[v.ne+1-i] ] = r[ i%r.ne+1 ];
    }

    return sparse( symmetric( m' * m ) );

};

ssi = function( A; n ) {

    # This function uses subspace iteration to compute the first
    # `n' eigenvalues and eigenvectors of `A'.

    local ( V; B; lambda; tmp; Afac; Vb; e; d; i; r );

    V = A[;1:n];
    B = diag( fill( A.nr; 1 ) );

    lambda = diag( transform(A;V) ) / diag( transform(B;V) );

    Afac = chol( A );
    while ( 1 ) {
        Vb = solve( Afac; V );
        e = eig( transform(A;Vb); transform(B;Vb) );
        d = e.values;
        tmp = abs( lambda - d )/ d;
        r = tmp[ imax( tmp ) ];
        lambda = e.values;
        V = Vb*e.vectors;
        if ( r < 1.0e-6 ) { break; };
    }

    return { values=e.values; vectors=V };

};

r = (
    0.22416506718106805,
    0.26479003963283732,
    0.76919488830920069,
    0.90072274622541049,
    0.089192990720827592,
    0.36610337829501527,
    0.27111135249543533,
    0.69340546694277105,
    0.32150303959916487,
    0.020990991509049661,
    0.60930041997195239,
    0.39150595217547657,
    0.049584821820997084,
    0.36070557793635205,
    0.047469744480899419,
    0.27127325594018831,
    0.3585841797099375,
    0.30496160001725031,
    0.63219573238501126,
    0.024773588881256799,
    0.065283286881299363,
    0.39865082194965834,
    0.83662774266518081,
    0.998176271560684,
    0.53930278194104453,
    0.90322313360088646,
    0.63726972259453951,
    0.050166096561665693,
    0.87081706378181323,
    0.039887622483022335,
    0.32569793394100754,
    0.094982130497220027,
    0.30467766211585962,
    0.094892822250208272,
    0.99570487718829181,
    0.39387065330234855,
    0.4609962005452235,
    0.26681622921806586,
    0.08727611977945833,
    0.78249924014438843,
    0.28780722072711551,
    0.69657653975141076,
    0.17400519185420368,
    0.33739204254811261,
    0.057282117222101483,
    0.2214749363351031,
    0.60866529895396215,
    0.41586629739770026,
    0.52643653635235343,
    0.24086103087331218,
    0.44063988627895706,
    0.5917198232336528,
    0.63951185282297052,
    0.27726762894413787,
    0.58989609432867551,
    0.17881463429835376,
    0.18049076207936313,
    0.22716581692321497,
    0.22898073132568073,
    0.051307825395515108,
    0.26705343940623733,
    0.5546786652666883,
    0.14628995635839642,
    0.57173110198775823,
    0.64957148751689653,
    0.14199483308102695,
    0.96560175529010672,
    0.1105676875964588,
    0.40881106229909281,
    0.052877874603903793,
    0.89306692820650846,
    0.69661828349186961,
    0.74945441482097586,
    0.067072119595050869,
    0.03401032557432089,
    0.80673653204307727,
    0.28854705639581524,
    0.64267562452828308,
    0.2226028289751163,
    0.81498359274816867,
    0.88353665540159521,
    0.6632427157197347,
    0.40670341598182147,
    0.52304850822456572,
    0.94051034466387251,
    0.99659951031049687,
    0.70186314252291948,
    0.12100110627757436,
    0.22376532676805058,
    0.93084387384860023,
    0.17230893213875076,
    0.49081876663994917,
    0.48552253911528853,
    0.31859888849714718,
    0.062549868162046121,
    0.13509402616652383,
    0.46059372157817413,
    0.028151622986491594,
    0.24566171422864391,
    0.86940478387726694
);

ssi( sparse_rand( 660; 0.01 ); 66 );