#include "Rocmop.h"
#include "roccom.h"
#include "Pane.h"
#include "Rocblas.h"
#include "Rocmap.h"
#include "geometry.h"
#include "Rocsurf.h"
#include "Generic_element_2.h"

Include dependency graph for smooth_medial.C:

Macros
#define	EXPORT_DEBUG_INFO 1

Functions
double	SQR (double x)

void	dsytrd3 (double A[3][3], double Q[3][3], double d[3], double e[2])

int	dsyevq3 (double A[3][3], double Q[3][3], double w[3])

Variables
static const double	half_pi = 1.5707963267949

static const double	pi = 3.14159265358979

Macro Definition Documentation

#define EXPORT_DEBUG_INFO 1

Definition at line 40 of file smooth_medial.C.

Function Documentation

int dsyevq3	(	double	A[3][3],
		double	Q[3][3],
		double	w[3]
	)

Definition at line 124 of file smooth_medial.C.

References dsytrd3(), i, k, s, SQR(), and sqrt().

Referenced by FaceOffset_3::compute_eigenvectors(), and Rocmop::compute_eigenvectors().

 {
   double e[3];                   // The third element is used only as temporary workspace
   double g, r, p, f, b, s, c, t; // Intermediate storage
   int nIter;
   int m;
 
   // Transform A to real tridiagonal form by the Householder method
   dsytrd3(A, Q, w, e);
   
   // Calculate eigensystem of the remaining real symmetric tridiagonal matrix
   // with the QL method
   //
   // Loop over all off-diagonal elements
   for (int l=0; l < 2; l++)
   {
     nIter = 0;
     for (;;)
     {
       // Check for convergence and exit iteration loop if off-diagonal
       // element e(l) is zero
       for (m=l; m < 2; m++)
       {
         g = fabs(w[m])+fabs(w[m+1]);
         if (fabs(e[m]) + g == g)
           break;
       }
       if (m == l)
         break;
       
       if (nIter++ >= 30)
         return -1;
 
       // Calculate g = d_m - k
       g = (w[l+1] - w[l]) / (e[l] + e[l]);
       r = sqrt(SQR(g) + 1.0);
       if (g > 0)
         g = w[m] - w[l] + e[l]/(g + r);
       else
         g = w[m] - w[l] + e[l]/(g - r);
 
       s = c = 1.0;
       p = 0.0;
       for (int i=m-1; i >= l; i--)
       {
         f = s * e[i];
         b = c * e[i];
         if (fabs(f) > fabs(g))
         {
           c      = g / f;
           r      = sqrt(SQR(c) + 1.0);
           e[i+1] = f * r;
           c     *= (s = 1.0/r);
         }
         else
         {
           s      = f / g;
           r      = sqrt(SQR(s) + 1.0);
           e[i+1] = g * r;
           s     *= (c = 1.0/r);
         }
         
         g = w[i+1] - p;
         r = (w[i] - g)*s + 2.0*c*b;
         p = s * r;
         w[i+1] = g + p;
         g = c*r - b;
 
         // Form eigenvectors
         for (int k=0; k < 3; k++)
         {
           t = Q[k][i+1];
           Q[k][i+1] = s*Q[k][i] + c*t;
           Q[k][i]   = c*Q[k][i] - s*t;
         }
       }
       w[l] -= p;
       e[l]  = g;
       e[m]  = 0.0;
     }
   }
 
   return 0;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void dsytrd3	(	double	A[3][3],
		double	Q[3][3],
		double	d[3],
		double	e[2]
	)

Definition at line 48 of file smooth_medial.C.

References i, j, q, SQR(), and sqrt().

Referenced by dsyevq3().

 {
   double u[3], q[3];
   double omega, f;
   double K, h, g;
   
   // Initialize Q to the identitity matrix
   for (int i=0; i < 3; i++)
   {
     Q[i][i] = 1.0;
     for (int j=0; j < i; j++)
       Q[i][j] = Q[j][i] = 0.0;
   }
 
   // Bring first row and column to the desired form 
   h = SQR(A[0][1]) + SQR(A[0][2]);
   if (A[0][1] > 0)
     g = -sqrt(h);
   else
     g = sqrt(h);
   e[0] = g;
   f    = g * A[0][1];
   u[1] = A[0][1] - g;
   u[2] = A[0][2];
   
   omega = h - f;
   if (omega > 0.0)
   {
     omega = 1.0 / omega;
     K     = 0.0;
     for (int i=1; i < 3; i++)
     {
       f    = A[1][i] * u[1] + A[i][2] * u[2];
       q[i] = omega * f;                  // p
       K   += u[i] * f;                   // u* A u
     }
     K *= 0.5 * SQR(omega);
 
     for (int i=1; i < 3; i++)
       q[i] = q[i] - K * u[i];
     
     d[0] = A[0][0];
     d[1] = A[1][1] - 2.0*q[1]*u[1];
     d[2] = A[2][2] - 2.0*q[2]*u[2];
     
     // Store inverse Householder transformation in Q
     for (int j=1; j < 3; j++)
     {
       f = omega * u[j];
       for (int i=1; i < 3; i++)
         Q[i][j] = Q[i][j] - f*u[i];
     }
 
     // Calculate updated A[1][2] and store it in e[1]
     e[1] = A[1][2] - q[1]*u[2] - u[1]*q[2];
   }
   else
   {
     for (int i=0; i < 3; i++)
       d[i] = A[i][i];
     e[1] = A[1][2];
   }
 }