src/eigsolv/SortingTools.cpp

Go to the documentation of this file.
00001 //**************************************************************************
00002 //
00003 //                                 NOTICE
00004 //
00005 // This software is a result of the research described in the report
00006 //
00007 // " A comparison of algorithms for modal analysis in the absence 
00008 //   of a sparse direct method", P. Arbenz, R. Lehoucq, and U. Hetmaniuk,
00009 //  Sandia National Laboratories, Technical report SAND2003-1028J.
00010 //
00011 // It is based on the Epetra, AztecOO, and ML packages defined in the Trilinos
00012 // framework ( http://software.sandia.gov/trilinos/ ).
00013 //
00014 // The distribution of this software follows also the rules defined in Trilinos.
00015 // This notice shall be marked on any reproduction of this software, in whole or
00016 // in part.
00017 //
00018 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00019 // license for use of this work by or on behalf of the U.S. Government.
00020 //
00021 // This program is distributed in the hope that it will be useful, but
00022 // WITHOUT ANY WARRANTY; without even the implied warranty of
00023 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00024 //
00025 // Code Authors: U. Hetmaniuk (ulhetma@sandia.gov), R. Lehoucq (rblehou@sandia.gov)
00026 //
00027 //**************************************************************************
00028 
00029 #include "SortingTools.h"
00030 
00031 
00032 int SortingTools::sortScalars(int n, double *y, int *perm) const {
00033 
00034     // Sort a vector into increasing order of algebraic values
00035     //
00036     // Input:
00037     //
00038     // n    (integer ) = Size of the array (input)
00039     // y    (double* ) = Array of length n to be sorted (input/output)
00040     // perm (integer*) = Array of length n with the permutation (input/output)
00041     //                   Optional argument
00042 
00043     int i, j;
00044     int igap = n / 2;
00045 
00046     if (igap == 0) {
00047         if ((n > 0) && (perm != 0)) {
00048             perm[0] = 0;
00049         }
00050         return 0;
00051     }
00052 
00053     if (perm) {
00054         for (i = 0; i < n; ++i)
00055             perm[i] = i;
00056     }
00057 
00058     while (igap > 0) {
00059         for (i=igap; i<n; ++i) {
00060             for (j=i-igap; j>=0; j-=igap) {
00061                 if (y[j] > y[j+igap]) {
00062                     double tmpD = y[j];
00063                     y[j] = y[j+igap];
00064                     y[j+igap] = tmpD;
00065                     if (perm) {
00066                         int tmpI = perm[j];
00067                         perm[j] = perm[j+igap];
00068                         perm[j+igap] = tmpI;
00069                     }
00070                 }
00071                 else {
00072                     break;
00073                 }
00074             }
00075         }
00076         igap = igap / 2;
00077     }
00078 
00079     return 0;
00080 
00081 }
00082 
00083 
00084 int SortingTools::sortScalars_Vectors(int num, double *lambda, double *Q, int ldQ) const {
00085 
00086     // This routines sorts the scalars (stored in lambda) in ascending order.
00087     // The associated vectors (stored in Q) are accordingly ordered.
00088     // One vector is of length ldQ.
00089     // Q must be of size ldQ * num.
00090 
00091     int info = 0;
00092     int i, j;
00093 
00094     int igap = num / 2;
00095 
00096     if ((Q) && (ldQ > 0)) {
00097         double *vec = new double[ldQ];
00098         double tmp;
00099         while (igap > 0) {
00100             for (i=igap; i < num; ++i) {
00101                 for (j=i-igap; j>=0; j-=igap) {
00102                     if (lambda[j] > lambda[j+igap]) {
00103                         tmp = lambda[j];
00104                         lambda[j] = lambda[j+igap];
00105                         lambda[j+igap] = tmp;
00107                         memcpy(vec, Q + j*ldQ, ldQ*sizeof(double));
00108                         memcpy(Q + j*ldQ, Q + (j+igap)*ldQ, ldQ*sizeof(double));
00109                         memcpy(Q + (j+igap)*ldQ, vec, ldQ*sizeof(double));
00110                     } 
00111                     else {
00112                         break;
00113                     }
00114                 }
00115             }
00116             igap = igap / 2;
00117         } // while (igap > 0)
00118         delete[] vec;
00119     } // if ((Q) && (ldQ > 0))
00120     else {
00121         while (igap > 0) {
00122             for (i=igap; i < num; ++i) {
00123                 for (j=i-igap; j>=0; j-=igap) {
00124                     if (lambda[j] > lambda[j+igap]) {
00125                         // Swap two scalars
00126                         double tmp = lambda[j];
00127                         lambda[j] = lambda[j+igap];
00128                         lambda[j+igap] = tmp;
00129                     } 
00130                     else {
00131                         break;
00132                     }
00133                 }
00134             }
00135             igap = igap / 2;
00136         } // while (igap > 0)
00137     } // if ((Q) && (ldQ > 0))
00138 
00139     return info;
00140 
00141 }
00142 
00143 

Generated on Fri Oct 26 13:35:11 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7