SVD decomposition.
519 if ( max_mpn > 0 && m+n > max_mpn )
521 error_msg =
"SVD_decomposition() error: m+n > " +
NOMAD::itos ( max_mpn );
525 double * rv1 =
new double[n];
533 int i , j , k , l = 0, its , jj , nm = 0;
534 double s , f , h , tmp , c , x , y , z , absf , absg , absh;
536 const int NITER = 30;
539 for (i=0; i < n; ++i)
542 for (j=0; j < n ; ++j)
547 for ( i = 0 ; i < n ; ++i )
554 for ( k = i ; k < m ; ++k )
555 scale += fabs ( M[k][i] );
558 for ( k = i ; k < m ; ++k )
561 s += M[k][i] * M[k][i];
564 g = ( f >= 0.0 ) ? -fabs(sqrt(s)) : fabs(sqrt(s));
567 for ( j = l ; j < n ; ++j )
569 for ( s = 0.0 , k = i ; k < m ; ++k )
570 s += M[k][i] * M[k][j];
572 for ( k = i ; k < m ; ++k )
573 M[k][j] += f * M[k][i];
575 for ( k = i ; k < m ; ++k )
581 if ( i < m && i != nm1 )
583 for ( k = l ; k < n ; ++k )
584 scale += fabs ( M[i][k] );
587 for ( k = l ; k < n ; ++k )
590 s += M[i][k] * M[i][k];
593 g = ( f >= 0.0 ) ? -fabs(sqrt(s)) : fabs(sqrt(s));
596 for ( k = l ; k < n ; ++k )
597 rv1[k] = M[i][k] / h;
598 for ( j = l ; j < m ; ++j )
600 for ( s=0.0,k=l ; k < n ; ++k )
601 s += M[j][k] * M[i][k];
602 for ( k=l ; k < n ; ++k )
603 M[j][k] += s * rv1[k];
605 for ( k = l ; k < n ; ++k )
609 tmp = fabs ( W[i] ) + fabs ( rv1[i] );
610 norm = ( norm > tmp ) ? norm : tmp;
614 for ( i = nm1 ; i >= 0 ; --i )
620 for ( j = l ; j < n ; ++j )
621 V[j][i] = ( M[i][j] / M[i][l] ) / g;
622 for ( j = l ; j < n ; ++j )
624 for ( s = 0.0 , k = l ; k < n ; ++k )
625 s += M[i][k] * V[k][j];
626 for ( k = l ; k < n ; ++k )
627 V[k][j] += s * V[k][i];
630 for ( j = l ; j < n ; ++j )
631 V[i][j] = V[j][i] = 0.0;
639 for ( i = ( ( m < n ) ? m : n ) - 1 ; i >= 0 ; --i )
643 for ( j = l ; j < n ; ++j )
648 for ( j = l ; j < n ; ++j )
650 for ( s = 0.0 , k = l ; k < m ; ++k )
651 s += M[k][i] * M[k][j];
652 f = ( s / M[i][i] ) * g;
653 for ( k = i ; k < m ; ++k )
654 M[k][j] += f * M[k][i];
656 for ( j = i ; j < m ; ++j )
660 for ( j = i ; j < m ; ++j )
666 for ( k = nm1 ; k >= 0 ; --k )
668 for ( its = 1 ; its <= NITER ; its++ )
671 for ( l = k ; l >= 0 ; l-- )
674 if ( nm < 0 || fabs ( rv1[l]) + norm == norm )
679 if ( fabs ( W[nm] ) + norm == norm )
686 for ( i = l ; i <= k ; i++ )
690 if ( fabs(f) + norm == norm )
696 h = ( absf > absg ) ?
697 absf * sqrt ( 1.0 + pow ( absg/absf , 2.0 ) ) :
698 ( ( absg==0 ) ? 0.0 : absg * sqrt ( 1.0 + pow ( absf/absg , 2.0 ) ) );
704 for ( j = 0 ; j < m ; ++j )
708 M[j][nm] = y * c + z * s;
709 M[j][ i] = z * c - y * s;
719 for ( j = 0 ; j < n ; j++ )
726 error_msg =
"SVD_decomposition() error: no convergence in " +
736 f = ( (y-z) * (y+z) + (g-h) * (g+h) ) / ( 2.0 * h * y );
740 absf * sqrt ( 1.0 + pow ( 1.0/absf , 2.0 ) ) :
741 sqrt ( 1.0 + pow ( absf , 2.0 ) );
743 f = ( (x-z) * (x+z) +
744 h * ( ( y / ( f + ( (f >= 0)? fabs(g) : -fabs(g) ) ) ) - h ) ) / x;
747 for ( j = l ; j <= nm ; ++j )
757 z = ( absf > absh ) ?
758 absf * sqrt ( 1.0 + pow ( absh/absf , 2.0 ) ) :
759 ( ( absh==0 ) ? 0.0 : absh * sqrt ( 1.0 + pow ( absf/absh , 2.0 ) ) );
768 for ( jj = 0 ; jj < n ; ++jj )
772 V[jj][j] = x * c + z * s;
773 V[jj][i] = z * c - x * s;
778 z = ( absf > absh ) ?
779 absf * sqrt ( 1.0 + pow ( absh/absf , 2.0 ) ) :
780 ( ( absh==0 ) ? 0.0 : absh * sqrt ( 1.0 + pow ( absf/absh , 2.0 ) ) );
792 for ( jj = 0 ; jj < m ; ++jj )
796 M[jj][j] = y * c + z * s;
797 M[jj][i] = z * c - y * s;
std::string itos(const int i)
Transform an integer into a string.