Mercurial > hg > octave-nkf

comparison liboctave/CMatrix.cc @ 4153:6b96ce9f5743

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-11-06 20:38:49 by jwe]
author jwe
date Wed, 06 Nov 2002 20:38:50 +0000
parents 7d9bda865012
children 5719210fff4c
comparison
equal deleted inserted replaced
4152:f14251d33b01 4153:6b96ce9f5743
984 const Complex *in (data ()); 984 const Complex *in (data ());
985 Complex *out (retval.fortran_vec ()); 985 Complex *out (retval.fortran_vec ());
986 986
987 for (size_t i = 0; i < nsamples; i++) 987 for (size_t i = 0; i < nsamples; i++)
988 { 988 {
989 OCTAVE_QUIT;
990
989 octave_fftw::fft (&in[npts * i], &out[npts * i], npts); 991 octave_fftw::fft (&in[npts * i], &out[npts * i], npts);
990 } 992 }
991 993
992 return retval; 994 return retval;
993 } 995 }
1016 const Complex *in (data ()); 1018 const Complex *in (data ());
1017 Complex *out (retval.fortran_vec ()); 1019 Complex *out (retval.fortran_vec ());
1018 1020
1019 for (size_t i = 0; i < nsamples; i++) 1021 for (size_t i = 0; i < nsamples; i++)
1020 { 1022 {
1023 OCTAVE_QUIT;
1024
1021 octave_fftw::ifft (&in[npts * i], &out[npts * i], npts); 1025 octave_fftw::ifft (&in[npts * i], &out[npts * i], npts);
1022 } 1026 }
1023 1027
1024 return retval; 1028 return retval;
1025 } 1029 }
1084 Complex *tmp_data = retval.fortran_vec (); 1088 Complex *tmp_data = retval.fortran_vec ();
1085 1089
1086 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1090 F77_FUNC (cffti, CFFTI) (npts, pwsave);
1087 1091
1088 for (int j = 0; j < nsamples; j++) 1092 for (int j = 0; j < nsamples; j++)
1089 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave); 1093 {
1094 OCTAVE_QUIT;
1095
1096 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
1097 }
1090 1098
1091 return retval; 1099 return retval;
1092 } 1100 }
1093 1101
1094 ComplexMatrix 1102 ComplexMatrix
1121 Complex *tmp_data = retval.fortran_vec (); 1129 Complex *tmp_data = retval.fortran_vec ();
1122 1130
1123 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1131 F77_FUNC (cffti, CFFTI) (npts, pwsave);
1124 1132
1125 for (int j = 0; j < nsamples; j++) 1133 for (int j = 0; j < nsamples; j++)
1126 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave); 1134 {
1135 OCTAVE_QUIT;
1136
1137 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
1138 }
1127 1139
1128 for (int j = 0; j < npts*nsamples; j++) 1140 for (int j = 0; j < npts*nsamples; j++)
1129 tmp_data[j] = tmp_data[j] / static_cast<double> (npts); 1141 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
1130 1142
1131 return retval; 1143 return retval;
1161 Complex *tmp_data = retval.fortran_vec (); 1173 Complex *tmp_data = retval.fortran_vec ();
1162 1174
1163 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1175 F77_FUNC (cffti, CFFTI) (npts, pwsave);
1164 1176
1165 for (int j = 0; j < nsamples; j++) 1177 for (int j = 0; j < nsamples; j++)
1166 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave); 1178 {
1179 OCTAVE_QUIT;
1180
1181 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
1182 }
1167 1183
1168 npts = nc; 1184 npts = nc;
1169 nsamples = nr; 1185 nsamples = nr;
1170 nn = 4*npts+15; 1186 nn = 4*npts+15;
1171 1187
1177 1193
1178 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1194 F77_FUNC (cffti, CFFTI) (npts, pwsave);
1179 1195
1180 for (int j = 0; j < nsamples; j++) 1196 for (int j = 0; j < nsamples; j++)
1181 { 1197 {
1198 OCTAVE_QUIT;
1199
1182 for (int i = 0; i < npts; i++) 1200 for (int i = 0; i < npts; i++)
1183 prow[i] = tmp_data[i*nr + j]; 1201 prow[i] = tmp_data[i*nr + j];
1184 1202
1185 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); 1203 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
1186 1204
1221 Complex *tmp_data = retval.fortran_vec (); 1239 Complex *tmp_data = retval.fortran_vec ();
1222 1240
1223 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1241 F77_FUNC (cffti, CFFTI) (npts, pwsave);
1224 1242
1225 for (int j = 0; j < nsamples; j++) 1243 for (int j = 0; j < nsamples; j++)
1226 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave); 1244 {
1245 OCTAVE_QUIT;
1246
1247 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
1248 }
1227 1249
1228 for (int j = 0; j < npts*nsamples; j++) 1250 for (int j = 0; j < npts*nsamples; j++)
1229 tmp_data[j] = tmp_data[j] / static_cast<double> (npts); 1251 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
1230 1252
1231 npts = nc; 1253 npts = nc;
1240 1262
1241 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1263 F77_FUNC (cffti, CFFTI) (npts, pwsave);
1242 1264
1243 for (int j = 0; j < nsamples; j++) 1265 for (int j = 0; j < nsamples; j++)
1244 { 1266 {
1267 OCTAVE_QUIT;
1268
1245 for (int i = 0; i < npts; i++) 1269 for (int i = 0; i < npts; i++)
1246 prow[i] = tmp_data[i*nr + j]; 1270 prow[i] = tmp_data[i*nr + j];
1247 1271
1248 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); 1272 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
1249 1273
1935 // inverse scaling (diagonal transformation) 1959 // inverse scaling (diagonal transformation)
1936 for (int i = 0; i < nc; i++) 1960 for (int i = 0; i < nc; i++)
1937 for (int j = 0; j < nc; j++) 1961 for (int j = 0; j < nc; j++)
1938 retval(i,j) *= dscale(i) / dscale(j); 1962 retval(i,j) *= dscale(i) / dscale(j);
1939 1963
1964 OCTAVE_QUIT;
1965
1940 // construct balancing permutation vector 1966 // construct balancing permutation vector
1941 Array<int> ipermute (nc); 1967 Array<int> ipermute (nc);
1942 for (int i = 0; i < nc; i++) 1968 for (int i = 0; i < nc; i++)
1943 ipermute(i) = i; // initialize to identity permutation 1969 ipermute(i) = i; // initialize to identity permutation
1944 1970
1962 1988
1963 // construct inverse balancing permutation vector 1989 // construct inverse balancing permutation vector
1964 Array<int> invpvec (nc); 1990 Array<int> invpvec (nc);
1965 for (int i = 0; i < nc; i++) 1991 for (int i = 0; i < nc; i++)
1966 invpvec(ipermute(i)) = i; // Thanks to R. A. Lippert for this method 1992 invpvec(ipermute(i)) = i; // Thanks to R. A. Lippert for this method
1993
1994 OCTAVE_QUIT;
1967 1995
1968 ComplexMatrix tmpMat = retval; 1996 ComplexMatrix tmpMat = retval;
1969 for (int i = 0; i < nc; i++) 1997 for (int i = 0; i < nc; i++)
1970 for (int j = 0; j < nc; j++) 1998 for (int j = 0; j < nc; j++)
1971 retval(i,j) = tmpMat(invpvec(i),invpvec(j)); 1999 retval(i,j) = tmpMat(invpvec(i),invpvec(j));