home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / cholp_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  2.4 KB  |  114 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Cholesky factorization with pivoting of a pos. 
  4. //              semi-definite matrix.
  5.  
  6. // Syntax:      CL = cholp ( A )
  7.  
  8. // Description:
  9.  
  10. //      Cholp returns R and a permutation matrix P such that 
  11.  
  12. //                      R'*R = P'*A*P
  13.  
  14. //      Only the upper triangular part of A is used. In addition the
  15. //      index I of the last positive diagonal element of R.  The first
  16. //      I rows of R contain the Cholesky factor of A.
  17.  
  18. //      This routine is based on the LINPACK routine CCHDC.  It works
  19. //      for both real and complex matrices. 
  20. //
  21. //      Reference:
  22. //      G.H. Golub and C.F. Van Loan, Matrix Computations, Second
  23. //      Edition, Johns Hopkins University Press, Baltimore, Maryland,
  24. //      1989, sec. 4.2.9.
  25.  
  26. //    This file is a translation of cholp.m from version 2.0 of
  27. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  28. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  29. //-------------------------------------------------------------------//
  30.  
  31. cholp = function ( A, pivot )
  32. {
  33.   local (A, pivot)
  34.  
  35.   if (exist (pivot)) { piv = pivot; else piv = 1; }
  36.  
  37.   n = A.nc;
  38.   pp = 1:n;
  39.   I = [];
  40.  
  41.   for (k in 1:n)
  42.   {
  43.     if (piv)
  44.     {
  45.       d = diag(A); 
  46.       big = max (d[k:n]);
  47.       m = maxi (d[k:n]);
  48.       m = m+k-1;
  49.     else 
  50.       big = A[k;k];  
  51.       m = k;
  52.     }
  53.     if (big <= 0) 
  54.     { 
  55.       I = k-1; 
  56.       break;
  57.     }
  58.  
  59.     // Symmetric row/column permutations.
  60.     if (m != k)
  61.     {
  62.       j = 1:k-1;
  63.       if (all(all(j)))
  64.       {
  65.         temp = A[j;k]; 
  66.         A[j;k] = A[j;m]; 
  67.         A[j;m] = temp;
  68.       }
  69.       temp = A[k;k]; 
  70.       A[k;k] = A[m;m]; 
  71.       A[m;m] = temp;
  72.       A[k;m] = conj(A[k;m]);
  73.       j = k+1:m-1;
  74.       if (all(all(j)))
  75.       {
  76.         temp = A[k;j]'; 
  77.         A[k;j] = A[j;m]'; 
  78.         A[j;m] = temp;
  79.       }
  80.       j = m+1:n;
  81.       if (all(all(j)))
  82.       {
  83.         temp = A[k;j]; 
  84.         A[k;j] = A[m;j]; 
  85.         A[m;j] = temp;
  86.       }
  87.       pp[k, m] = pp[m, k];
  88.     }
  89.  
  90.     A[k;k] = sqrt( A[k;k] );
  91.     if (k == n) { break }
  92.     A[k; k+1:n] = A[k; k+1:n] / A[k;k];
  93.  
  94.     // For simplicity update the whole of the remaining submatrix
  95.     // (rather than just the upper triangle).
  96.  
  97.     j = k+1:n;
  98.     A[j;j] = A[j;j] - A[k;j]'*A[k;j];
  99.  
  100.   }  // Main k loop
  101.  
  102.   R = triu(A);
  103.   if (isempty(I)) { I = n; }
  104.   if (!piv)
  105.   {
  106.     P = I;
  107.   else
  108.    P = eye(n,n); 
  109.    P = P[;pp];
  110.   }
  111.  
  112.   return << R = R; P = P; I = I >>;
  113. };
  114.