ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/math/SquareMatrix.hpp
(Generate patch)

Comparing trunk/src/math/SquareMatrix.hpp (file contents):
Revision 137 by tim, Thu Oct 21 21:31:39 2004 UTC vs.
Revision 151 by tim, Mon Oct 25 22:46:19 2004 UTC

# Line 55 | Line 55 | namespace oopse {
55                          data_[i][j] = 0.0;
56               }
57  
58 +            /** Constructs and initializes every element of this matrix to a scalar */
59 +            SquareMatrix(Real s) : RectMatrix<Real, Dim, Dim>(s){
60 +            }
61 +
62 +            /** Constructs and initializes from an array */
63 +            SquareMatrix(Real* array) : RectMatrix<Real, Dim, Dim>(array){
64 +            }
65 +
66 +
67              /** copy constructor */
68              SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
69              }
# Line 199 | Line 208 | namespace oopse {
208      // normalized.
209      template<typename Real, int Dim>
210      int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
211 <                                  SquareMatrix<Real, Dim>& v) {
212 <      const int n = Dim;  
213 <      int i, j, k, iq, ip, numPos;
214 <      Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
215 <      Real bspace[4], zspace[4];
216 <      Real *b = bspace;
217 <      Real *z = zspace;
211 >                                        SquareMatrix<Real, Dim>& v) {
212 >        const int n = Dim;  
213 >        int i, j, k, iq, ip, numPos;
214 >        Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
215 >        Real bspace[4], zspace[4];
216 >        Real *b = bspace;
217 >        Real *z = zspace;
218  
219 <      // only allocate memory if the matrix is large
220 <      if (n > 4)
221 <        {
222 <        b = new Real[n];
214 <        z = new Real[n];
219 >        // only allocate memory if the matrix is large
220 >        if (n > 4) {
221 >            b = new Real[n];
222 >            z = new Real[n];
223          }
224  
225 <      // initialize
226 <      for (ip=0; ip<n; ip++)
227 <        {
228 <        for (iq=0; iq<n; iq++)
229 <          {
230 <          v(ip, iq) = 0.0;
223 <          }
224 <        v(ip, ip) = 1.0;
225 >        // initialize
226 >        for (ip=0; ip<n; ip++) {
227 >            for (iq=0; iq<n; iq++) {
228 >                v(ip, iq) = 0.0;
229 >            }
230 >            v(ip, ip) = 1.0;
231          }
232 <      for (ip=0; ip<n; ip++)
233 <        {
234 <        b[ip] = w[ip] = a(ip, ip);
229 <        z[ip] = 0.0;
232 >        for (ip=0; ip<n; ip++) {
233 >            b[ip] = w[ip] = a(ip, ip);
234 >            z[ip] = 0.0;
235          }
236  
237 <      // begin rotation sequence
238 <      for (i=0; i<VTK_MAX_ROTATIONS; i++)
239 <        {
240 <        sm = 0.0;
241 <        for (ip=0; ip<n-1; ip++)
242 <          {
243 <          for (iq=ip+1; iq<n; iq++)
239 <            {
240 <            sm += fabs(a(ip, iq));
237 >        // begin rotation sequence
238 >        for (i=0; i<VTK_MAX_ROTATIONS; i++) {
239 >            sm = 0.0;
240 >            for (ip=0; ip<n-1; ip++) {
241 >                for (iq=ip+1; iq<n; iq++) {
242 >                    sm += fabs(a(ip, iq));
243 >                }
244              }
245 <          }
246 <        if (sm == 0.0)
247 <          {
245 <          break;
246 <          }
245 >            if (sm == 0.0) {
246 >                break;
247 >            }
248  
249 <        if (i < 3)                                // first 3 sweeps
250 <          {
251 <          tresh = 0.2*sm/(n*n);
252 <          }
253 <        else
253 <          {
254 <          tresh = 0.0;
255 <          }
249 >            if (i < 3) {                                // first 3 sweeps
250 >                tresh = 0.2*sm/(n*n);
251 >            } else {
252 >                tresh = 0.0;
253 >            }
254  
255 <        for (ip=0; ip<n-1; ip++)
256 <          {
257 <          for (iq=ip+1; iq<n; iq++)
260 <            {
261 <            g = 100.0*fabs(a(ip, iq));
255 >            for (ip=0; ip<n-1; ip++) {
256 >                for (iq=ip+1; iq<n; iq++) {
257 >                    g = 100.0*fabs(a(ip, iq));
258  
259 <            // after 4 sweeps
260 <            if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
261 <            && (fabs(w[iq])+g) == fabs(w[iq]))
262 <              {
263 <              a(ip, iq) = 0.0;
264 <              }
265 <            else if (fabs(a(ip, iq)) > tresh)
266 <              {
267 <              h = w[iq] - w[ip];
268 <              if ( (fabs(h)+g) == fabs(h))
269 <                {
270 <                t = (a(ip, iq)) / h;
271 <                }
272 <              else
273 <                {
274 <                theta = 0.5*h / (a(ip, iq));
275 <                t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
276 <                if (theta < 0.0)
277 <                  {
278 <                  t = -t;
279 <                  }
280 <                }
281 <              c = 1.0 / sqrt(1+t*t);
282 <              s = t*c;
283 <              tau = s/(1.0+c);
284 <              h = t*a(ip, iq);
285 <              z[ip] -= h;
286 <              z[iq] += h;
287 <              w[ip] -= h;
288 <              w[iq] += h;
289 <              a(ip, iq)=0.0;
290 <
291 <              // ip already shifted left by 1 unit
292 <              for (j = 0;j <= ip-1;j++)
293 <                {
294 <                VTK_ROTATE(a,j,ip,j,iq);
295 <                }
296 <              // ip and iq already shifted left by 1 unit
297 <              for (j = ip+1;j <= iq-1;j++)
298 <                {
299 <                VTK_ROTATE(a,ip,j,j,iq);
304 <                }
305 <              // iq already shifted left by 1 unit
306 <              for (j=iq+1; j<n; j++)
307 <                {
308 <                VTK_ROTATE(a,ip,j,iq,j);
309 <                }
310 <              for (j=0; j<n; j++)
311 <                {
312 <                VTK_ROTATE(v,j,ip,j,iq);
259 >                    // after 4 sweeps
260 >                    if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
261 >                        && (fabs(w[iq])+g) == fabs(w[iq])) {
262 >                        a(ip, iq) = 0.0;
263 >                    } else if (fabs(a(ip, iq)) > tresh) {
264 >                        h = w[iq] - w[ip];
265 >                        if ( (fabs(h)+g) == fabs(h)) {
266 >                            t = (a(ip, iq)) / h;
267 >                        } else {
268 >                            theta = 0.5*h / (a(ip, iq));
269 >                            t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
270 >                            if (theta < 0.0) {
271 >                                t = -t;
272 >                            }
273 >                        }
274 >                        c = 1.0 / sqrt(1+t*t);
275 >                        s = t*c;
276 >                        tau = s/(1.0+c);
277 >                        h = t*a(ip, iq);
278 >                        z[ip] -= h;
279 >                        z[iq] += h;
280 >                        w[ip] -= h;
281 >                        w[iq] += h;
282 >                        a(ip, iq)=0.0;
283 >
284 >                        // ip already shifted left by 1 unit
285 >                        for (j = 0;j <= ip-1;j++) {
286 >                            VTK_ROTATE(a,j,ip,j,iq);
287 >                        }
288 >                        // ip and iq already shifted left by 1 unit
289 >                        for (j = ip+1;j <= iq-1;j++) {
290 >                            VTK_ROTATE(a,ip,j,j,iq);
291 >                        }
292 >                        // iq already shifted left by 1 unit
293 >                        for (j=iq+1; j<n; j++) {
294 >                            VTK_ROTATE(a,ip,j,iq,j);
295 >                        }
296 >                        for (j=0; j<n; j++) {
297 >                            VTK_ROTATE(v,j,ip,j,iq);
298 >                        }
299 >                    }
300                  }
314 –              }
301              }
316 –          }
302  
303 <        for (ip=0; ip<n; ip++)
304 <          {
305 <          b[ip] += z[ip];
306 <          w[ip] = b[ip];
307 <          z[ip] = 0.0;
323 <          }
303 >            for (ip=0; ip<n; ip++) {
304 >                b[ip] += z[ip];
305 >                w[ip] = b[ip];
306 >                z[ip] = 0.0;
307 >            }
308          }
309  
310 <      //// this is NEVER called
311 <      if ( i >= VTK_MAX_ROTATIONS )
312 <        {
313 <           std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
330 <           return 0;
310 >        //// this is NEVER called
311 >        if ( i >= VTK_MAX_ROTATIONS ) {
312 >            std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
313 >            return 0;
314          }
315  
316 <      // sort eigenfunctions                 these changes do not affect accuracy
317 <      for (j=0; j<n-1; j++)                  // boundary incorrect
318 <        {
336 <        k = j;
337 <        tmp = w[k];
338 <        for (i=j+1; i<n; i++)                // boundary incorrect, shifted already
339 <          {
340 <          if (w[i] >= tmp)                   // why exchage if same?
341 <            {
342 <            k = i;
316 >        // sort eigenfunctions                 these changes do not affect accuracy
317 >        for (j=0; j<n-1; j++) {                  // boundary incorrect
318 >            k = j;
319              tmp = w[k];
320 +            for (i=j+1; i<n; i++) {                // boundary incorrect, shifted already
321 +                if (w[i] >= tmp) {                   // why exchage if same?
322 +                    k = i;
323 +                    tmp = w[k];
324 +                }
325              }
326 <          }
327 <        if (k != j)
328 <          {
329 <          w[k] = w[j];
330 <          w[j] = tmp;
331 <          for (i=0; i<n; i++)
332 <            {
333 <            tmp = v(i, j);
353 <            v(i, j) = v(i, k);
354 <            v(i, k) = tmp;
326 >            if (k != j) {
327 >                w[k] = w[j];
328 >                w[j] = tmp;
329 >                for (i=0; i<n; i++) {
330 >                    tmp = v(i, j);
331 >                    v(i, j) = v(i, k);
332 >                    v(i, k) = tmp;
333 >                }
334              }
356 –          }
335          }
336 <      // insure eigenvector consistency (i.e., Jacobi can compute vectors that
337 <      // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
338 <      // reek havoc in hyperstreamline/other stuff. We will select the most
339 <      // positive eigenvector.
340 <      int ceil_half_n = (n >> 1) + (n & 1);
341 <      for (j=0; j<n; j++)
342 <        {
343 <        for (numPos=0, i=0; i<n; i++)
344 <          {
345 <          if ( v(i, j) >= 0.0 )
368 <            {
369 <            numPos++;
336 >        // insure eigenvector consistency (i.e., Jacobi can compute vectors that
337 >        // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
338 >        // reek havoc in hyperstreamline/other stuff. We will select the most
339 >        // positive eigenvector.
340 >        int ceil_half_n = (n >> 1) + (n & 1);
341 >        for (j=0; j<n; j++) {
342 >            for (numPos=0, i=0; i<n; i++) {
343 >                if ( v(i, j) >= 0.0 ) {
344 >                    numPos++;
345 >                }
346              }
347 <          }
348 <    //    if ( numPos < ceil(double(n)/double(2.0)) )
349 <        if ( numPos < ceil_half_n)
350 <          {
351 <          for(i=0; i<n; i++)
376 <            {
377 <            v(i, j) *= -1.0;
347 >            //    if ( numPos < ceil(double(n)/double(2.0)) )
348 >            if ( numPos < ceil_half_n) {
349 >                for (i=0; i<n; i++) {
350 >                    v(i, j) *= -1.0;
351 >                }
352              }
379 –          }
353          }
354  
355 <      if (n > 4)
356 <        {
357 <        delete [] b;
385 <        delete [] z;
355 >        if (n > 4) {
356 >            delete [] b;
357 >            delete [] z;
358          }
359 <      return 1;
359 >        return 1;
360      }
361  
362  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines