Mercurial > octave-nkf

comparison liboctave/CMatrix.cc @ 1948:d7dec93d4b87

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-02-14 03:45:49 by jwe]
author jwe
date Wed, 14 Feb 1996 03:47:59 +0000
parents 1281a23a34dd
children ab6abe89aaa1
comparison
equal deleted inserted replaced
1947:5ab3c25a3cf9 1948:d7dec93d4b87
799 } 799 }
800 800
801 ComplexMatrix 801 ComplexMatrix
802 ComplexMatrix::inverse (int& info, double& rcond, int force) const 802 ComplexMatrix::inverse (int& info, double& rcond, int force) const
803 { 803 {
804 int nr = rows (); 804 ComplexMatrix retval;
805 int nc = cols (); 805
806 int len = length (); 806 int nr = rows ();
807 int nc = cols ();
808
807 if (nr != nc) 809 if (nr != nc)
808 { 810 (*current_liboctave_error_handler) ("inverse requires square matrix");
809 (*current_liboctave_error_handler) ("inverse requires square matrix");
810 return ComplexMatrix ();
811 }
812
813 info = 0;
814
815 int *ipvt = new int [nr];
816 Complex *z = new Complex [nr];
817 Complex *tmp_data = dup (data (), len);
818
819 F77_FCN (zgeco, ZGECO) (tmp_data, nr, nc, ipvt, rcond, z);
820
821 volatile double rcond_plus_one = rcond + 1.0;
822
823 if (rcond_plus_one == 1.0)
824 info = -1;
825
826 if (info == -1 && ! force)
827 {
828 copy (tmp_data, data (), len); // Restore contents.
829 }
830 else 811 else
831 { 812 {
832 Complex *dummy = 0; 813 info = 0;
833 814
834 F77_FCN (zgedi, ZGEDI) (tmp_data, nr, nc, ipvt, dummy, z, 1); 815 Array<int> ipvt (nr);
835 } 816 int *pipvt = ipvt.fortran_vec ();
836 817
837 delete [] ipvt; 818 Array<Complex> z (nr);
838 delete [] z; 819 Complex *pz = z.fortran_vec ();
839 820
840 return ComplexMatrix (tmp_data, nr, nc); 821 retval = *this;
822 Complex *tmp_data = retval.fortran_vec ();
823
824 F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nc, pipvt, rcond, pz));
825
826 if (f77_exception_encountered)
827 (*current_liboctave_error_handler) ("unrecoverable error in zgeco");
828 else
829 {
830 volatile double rcond_plus_one = rcond + 1.0;
831
832 if (rcond_plus_one == 1.0)
833 info = -1;
834
835 if (info == -1 && ! force)
836 retval = *this; // Restore contents.
837 else
838 {
839 Complex *dummy = 0;
840
841 F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nc, pipvt, dummy,
842 pz, 1));
843
844 if (f77_exception_encountered)
845 (*current_liboctave_error_handler)
846 ("unrecoverable error in zgedi");
847 }
848 }
849 }
850
851 return retval;
841 } 852 }
842 853
843 ComplexMatrix 854 ComplexMatrix
844 ComplexMatrix::pseudo_inverse (double tol) 855 ComplexMatrix::pseudo_inverse (double tol)
845 { 856 {
882 } 893 }
883 894
884 ComplexMatrix 895 ComplexMatrix
885 ComplexMatrix::fourier (void) const 896 ComplexMatrix::fourier (void) const
886 { 897 {
887 int nr = rows (); 898 ComplexMatrix retval;
888 int nc = cols (); 899
900 int nr = rows ();
901 int nc = cols ();
902
889 int npts, nsamples; 903 int npts, nsamples;
904
890 if (nr == 1 || nc == 1) 905 if (nr == 1 || nc == 1)
891 { 906 {
892 npts = nr > nc ? nr : nc; 907 npts = nr > nc ? nr : nc;
893 nsamples = 1; 908 nsamples = 1;
894 } 909 }
897 npts = nr; 912 npts = nr;
898 nsamples = nc; 913 nsamples = nc;
899 } 914 }
900 915
901 int nn = 4*npts+15; 916 int nn = 4*npts+15;
902 Complex *wsave = new Complex [nn]; 917
903 Complex *tmp_data = dup (data (), length ()); 918 Array<Complex> wsave (nn);
904 919 Complex *pwsave = wsave.fortran_vec ();
905 F77_FCN (cffti, CFFTI) (npts, wsave); 920
921 retval = *this;
922 Complex *tmp_data = retval.fortran_vec ();
923
924 F77_FCN (cffti, CFFTI) (npts, pwsave);
906 925
907 for (int j = 0; j < nsamples; j++) 926 for (int j = 0; j < nsamples; j++)
908 F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], wsave); 927 F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
909 928
910 delete [] wsave; 929 return retval;
911
912 return ComplexMatrix (tmp_data, nr, nc);
913 } 930 }
914 931
915 ComplexMatrix 932 ComplexMatrix
916 ComplexMatrix::ifourier (void) const 933 ComplexMatrix::ifourier (void) const
917 { 934 {
918 int nr = rows (); 935 ComplexMatrix retval;
919 int nc = cols (); 936
937 int nr = rows ();
938 int nc = cols ();
939
920 int npts, nsamples; 940 int npts, nsamples;
941
921 if (nr == 1 || nc == 1) 942 if (nr == 1 || nc == 1)
922 { 943 {
923 npts = nr > nc ? nr : nc; 944 npts = nr > nc ? nr : nc;
924 nsamples = 1; 945 nsamples = 1;
925 } 946 }
928 npts = nr; 949 npts = nr;
929 nsamples = nc; 950 nsamples = nc;
930 } 951 }
931 952
932 int nn = 4*npts+15; 953 int nn = 4*npts+15;
933 Complex *wsave = new Complex [nn]; 954
934 Complex *tmp_data = dup (data (), length ()); 955 Array<Complex> wsave (nn);
935 956 Complex *pwsave = wsave.fortran_vec ();
936 F77_FCN (cffti, CFFTI) (npts, wsave); 957
958 retval = *this;
959 Complex *tmp_data = retval.fortran_vec ();
960
961 F77_FCN (cffti, CFFTI) (npts, pwsave);
937 962
938 for (int j = 0; j < nsamples; j++) 963 for (int j = 0; j < nsamples; j++)
939 F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], wsave); 964 F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
940 965
941 for (int j = 0; j < npts*nsamples; j++) 966 for (int j = 0; j < npts*nsamples; j++)
942 tmp_data[j] = tmp_data[j] / (double) npts; 967 tmp_data[j] = tmp_data[j] / (double) npts;
943 968
944 delete [] wsave; 969 return retval;
945
946 return ComplexMatrix (tmp_data, nr, nc);
947 } 970 }
948 971
949 ComplexMatrix 972 ComplexMatrix
950 ComplexMatrix::fourier2d (void) const 973 ComplexMatrix::fourier2d (void) const
951 { 974 {
952 int nr = rows (); 975 ComplexMatrix retval;
953 int nc = cols (); 976
977 int nr = rows ();
978 int nc = cols ();
979
954 int npts, nsamples; 980 int npts, nsamples;
981
955 if (nr == 1 || nc == 1) 982 if (nr == 1 || nc == 1)
956 { 983 {
957 npts = nr > nc ? nr : nc; 984 npts = nr > nc ? nr : nc;
958 nsamples = 1; 985 nsamples = 1;
959 } 986 }
962 npts = nr; 989 npts = nr;
963 nsamples = nc; 990 nsamples = nc;
964 } 991 }
965 992
966 int nn = 4*npts+15; 993 int nn = 4*npts+15;
967 Complex *wsave = new Complex [nn]; 994
968 Complex *tmp_data = dup (data (), length ()); 995 Array<Complex> wsave (nn);
969 996 Complex *pwsave = wsave.fortran_vec ();
970 F77_FCN (cffti, CFFTI) (npts, wsave); 997
998 retval = *this;
999 Complex *tmp_data = retval.fortran_vec ();
1000
1001 F77_FCN (cffti, CFFTI) (npts, pwsave);
971 1002
972 for (int j = 0; j < nsamples; j++) 1003 for (int j = 0; j < nsamples; j++)
973 F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], wsave); 1004 F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
974
975 delete [] wsave;
976 1005
977 npts = nc; 1006 npts = nc;
978 nsamples = nr; 1007 nsamples = nr;
979 nn = 4*npts+15; 1008 nn = 4*npts+15;
980 wsave = new Complex [nn]; 1009
981 Complex *row = new Complex[npts]; 1010 wsave.resize (nn);
982 1011 pwsave = wsave.fortran_vec ();
983 F77_FCN (cffti, CFFTI) (npts, wsave); 1012
1013 Array<Complex> row (npts);
1014 Complex *prow = row.fortran_vec ();
1015
1016 F77_FCN (cffti, CFFTI) (npts, pwsave);
984 1017
985 for (int j = 0; j < nsamples; j++) 1018 for (int j = 0; j < nsamples; j++)
986 { 1019 {
987 for (int i = 0; i < npts; i++) 1020 for (int i = 0; i < npts; i++)
988 row[i] = tmp_data[i*nr + j]; 1021 prow[i] = tmp_data[i*nr + j];
989 1022
990 F77_FCN (cfftf, CFFTF) (npts, row, wsave); 1023 F77_FCN (cfftf, CFFTF) (npts, prow, pwsave);
991 1024
992 for (int i = 0; i < npts; i++) 1025 for (int i = 0; i < npts; i++)
993 tmp_data[i*nr + j] = row[i]; 1026 tmp_data[i*nr + j] = prow[i];
994 } 1027 }
995 1028
996 delete [] wsave; 1029 return retval;
997 delete [] row;
998
999 return ComplexMatrix (tmp_data, nr, nc);
1000 } 1030 }
1001 1031
1002 ComplexMatrix 1032 ComplexMatrix
1003 ComplexMatrix::ifourier2d (void) const 1033 ComplexMatrix::ifourier2d (void) const
1004 { 1034 {
1005 int nr = rows (); 1035 ComplexMatrix retval;
1006 int nc = cols (); 1036
1037 int nr = rows ();
1038 int nc = cols ();
1039
1007 int npts, nsamples; 1040 int npts, nsamples;
1041
1008 if (nr == 1 || nc == 1) 1042 if (nr == 1 || nc == 1)
1009 { 1043 {
1010 npts = nr > nc ? nr : nc; 1044 npts = nr > nc ? nr : nc;
1011 nsamples = 1; 1045 nsamples = 1;
1012 } 1046 }
1015 npts = nr; 1049 npts = nr;
1016 nsamples = nc; 1050 nsamples = nc;
1017 } 1051 }
1018 1052
1019 int nn = 4*npts+15; 1053 int nn = 4*npts+15;
1020 Complex *wsave = new Complex [nn]; 1054
1021 Complex *tmp_data = dup (data (), length ()); 1055 Array<Complex> wsave (nn);
1022 1056 Complex *pwsave = wsave.fortran_vec ();
1023 F77_FCN (cffti, CFFTI) (npts, wsave); 1057
1058 retval = *this;
1059 Complex *tmp_data = retval.fortran_vec ();
1060
1061 F77_FCN (cffti, CFFTI) (npts, pwsave);
1024 1062
1025 for (int j = 0; j < nsamples; j++) 1063 for (int j = 0; j < nsamples; j++)
1026 F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], wsave); 1064 F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
1027
1028 delete [] wsave;
1029 1065
1030 for (int j = 0; j < npts*nsamples; j++) 1066 for (int j = 0; j < npts*nsamples; j++)
1031 tmp_data[j] = tmp_data[j] / (double) npts; 1067 tmp_data[j] = tmp_data[j] / (double) npts;
1032 1068
1033 npts = nc; 1069 npts = nc;
1034 nsamples = nr; 1070 nsamples = nr;
1035 nn = 4*npts+15; 1071 nn = 4*npts+15;
1036 wsave = new Complex [nn]; 1072
1037 Complex *row = new Complex[npts]; 1073 wsave.resize (nn);
1038 1074 pwsave = wsave.fortran_vec ();
1039 F77_FCN (cffti, CFFTI) (npts, wsave); 1075
1076 Array<Complex> row (npts);
1077 Complex *prow = row.fortran_vec ();
1078
1079 F77_FCN (cffti, CFFTI) (npts, pwsave);
1040 1080
1041 for (int j = 0; j < nsamples; j++) 1081 for (int j = 0; j < nsamples; j++)
1042 { 1082 {
1043 for (int i = 0; i < npts; i++) 1083 for (int i = 0; i < npts; i++)
1044 row[i] = tmp_data[i*nr + j]; 1084 prow[i] = tmp_data[i*nr + j];
1045 1085
1046 F77_FCN (cfftb, CFFTB) (npts, row, wsave); 1086 F77_FCN (cfftb, CFFTB) (npts, prow, pwsave);
1047 1087
1048 for (int i = 0; i < npts; i++) 1088 for (int i = 0; i < npts; i++)
1049 tmp_data[i*nr + j] = row[i] / (double) npts; 1089 tmp_data[i*nr + j] = prow[i] / (double) npts;
1050 } 1090 }
1051 1091
1052 delete [] wsave; 1092 return retval;
1053 delete [] row;
1054
1055 return ComplexMatrix (tmp_data, nr, nc);
1056 } 1093 }
1057 1094
1058 ComplexDET 1095 ComplexDET
1059 ComplexMatrix::determinant (void) const 1096 ComplexMatrix::determinant (void) const
1060 { 1097 {
1086 retval = ComplexDET (d); 1123 retval = ComplexDET (d);
1087 } 1124 }
1088 else 1125 else
1089 { 1126 {
1090 info = 0; 1127 info = 0;
1091 int *ipvt = new int [nr]; 1128
1092 1129 Array<int> ipvt (nr);
1093 Complex *z = new Complex [nr]; 1130 int *pipvt = ipvt.fortran_vec ();
1094 Complex *tmp_data = dup (data (), length ()); 1131
1095 1132 Array<Complex> z (nr);
1096 F77_FCN (zgeco, ZGECO) (tmp_data, nr, nr, ipvt, rcond, z); 1133 Complex *pz = z.fortran_vec ();
1097 1134
1098 volatile double rcond_plus_one = rcond + 1.0; 1135 ComplexMatrix atmp = *this;
1099 if (rcond_plus_one == 1.0) 1136 Complex *tmp_data = atmp.fortran_vec ();
1100 { 1137
1101 info = -1; 1138 F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
1102 retval = ComplexDET (); 1139
1103 } 1140 if (f77_exception_encountered)
1141 (*current_liboctave_error_handler) ("unrecoverable error in zgeco");
1104 else 1142 else
1105 { 1143 {
1106 Complex d[2]; 1144 volatile double rcond_plus_one = rcond + 1.0;
1107 F77_FCN (zgedi, ZGEDI) (tmp_data, nr, nr, ipvt, d, z, 10); 1145
1108 retval = ComplexDET (d); 1146 if (rcond_plus_one == 1.0)
1109 } 1147 {
1110 1148 info = -1;
1111 delete [] tmp_data; 1149 retval = ComplexDET ();
1112 delete [] ipvt; 1150 }
1113 delete [] z; 1151 else
1152 {
1153 Complex d[2];
1154
1155 F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nr, pipvt, d, pz, 10));
1156
1157 if (f77_exception_encountered)
1158 (*current_liboctave_error_handler)
1159 ("unrecoverable error in dgedi");
1160 else
1161 retval = ComplexDET (d);
1162 }
1163 }
1114 } 1164 }
1115 1165
1116 return retval; 1166 return retval;
1117 } 1167 }
1118 1168
1157 { 1207 {
1158 ComplexMatrix retval; 1208 ComplexMatrix retval;
1159 1209
1160 int nr = rows (); 1210 int nr = rows ();
1161 int nc = cols (); 1211 int nc = cols ();
1162 int b_nr = b.rows (); 1212
1163 int b_nc = b.cols (); 1213 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
1164 if (nr == 0 || nc == 0 || nr != nc || nr != b_nr) 1214 (*current_liboctave_error_handler)
1165 { 1215 ("matrix dimension mismatch in solution of linear equations");
1166 (*current_liboctave_error_handler)
1167 ("matrix dimension mismatch in solution of linear equations");
1168 return ComplexMatrix ();
1169 }
1170
1171 info = 0;
1172 int *ipvt = new int [nr];
1173
1174 Complex *z = new Complex [nr];
1175 Complex *tmp_data = dup (data (), length ());
1176
1177 F77_FCN (zgeco, ZGECO) (tmp_data, nr, nr, ipvt, rcond, z);
1178
1179 volatile double rcond_plus_one = rcond + 1.0;
1180 if (rcond_plus_one == 1.0)
1181 {
1182 info = -2;
1183 }
1184 else 1216 else
1185 { 1217 {
1186 Complex *result = dup (b.data (), b.length ()); 1218 info = 0;
1187 1219
1188 for (int j = 0; j < b_nc; j++) 1220 Array<int> ipvt (nr);
1189 F77_FCN (zgesl, ZGESL) (tmp_data, nr, nr, ipvt, &result[nr*j], 0); 1221 int *pipvt = ipvt.fortran_vec ();
1190 1222
1191 retval = ComplexMatrix (result, b_nr, b_nc); 1223 Array<Complex> z (nr);
1192 } 1224 Complex *pz = z.fortran_vec ();
1193 1225
1194 delete [] tmp_data; 1226 ComplexMatrix atmp = *this;
1195 delete [] ipvt; 1227 Complex *tmp_data = atmp.fortran_vec ();
1196 delete [] z; 1228
1229 F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
1230
1231 if (f77_exception_encountered)
1232 (*current_liboctave_error_handler) ("unrecoverable error in zgeco");
1233 else
1234 {
1235 volatile double rcond_plus_one = rcond + 1.0;
1236
1237 if (rcond_plus_one == 1.0)
1238 {
1239 info = -2;
1240 }
1241 else
1242 {
1243 retval = b;
1244 Complex *result = retval.fortran_vec ();
1245
1246 int b_nc = b.cols ();
1247
1248 for (volatile int j = 0; j < b_nc; j++)
1249 {
1250 F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt,
1251 &result[nr*j], 0));
1252
1253 if (f77_exception_encountered)
1254 {
1255 (*current_liboctave_error_handler)
1256 ("unrecoverable error in dgesl");
1257
1258 break;
1259 }
1260 }
1261 }
1262 }
1263 }
1197 1264
1198 return retval; 1265 return retval;
1199 } 1266 }
1200 1267
1201 ComplexColumnVector 1268 ComplexColumnVector
1219 { 1286 {
1220 ComplexColumnVector retval; 1287 ComplexColumnVector retval;
1221 1288
1222 int nr = rows (); 1289 int nr = rows ();
1223 int nc = cols (); 1290 int nc = cols ();
1224 int b_len = b.length (); 1291
1225 if (nr == 0 || nc == 0 || nr != nc || nr != b_len) 1292 if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
1226 { 1293 (*current_liboctave_error_handler)
1227 (*current_liboctave_error_handler) 1294 ("matrix dimension mismatch in solution of linear equations");
1228 ("matrix dimension mismatch in solution of linear equations");
1229 return ComplexColumnVector ();
1230 }
1231
1232 info = 0;
1233 int *ipvt = new int [nr];
1234
1235 Complex *z = new Complex [nr];
1236 Complex *tmp_data = dup (data (), length ());
1237
1238 F77_FCN (zgeco, ZGECO) (tmp_data, nr, nr, ipvt, rcond, z);
1239
1240 volatile double rcond_plus_one = rcond + 1.0;
1241 if (rcond_plus_one == 1.0)
1242 {
1243 info = -2;
1244 }
1245 else 1295 else
1246 { 1296 {
1247 Complex *result = dup (b.data (), b_len); 1297 info = 0;
1248 1298
1249 F77_FCN (zgesl, ZGESL) (tmp_data, nr, nr, ipvt, result, 0); 1299 Array<int> ipvt (nr);
1250 1300 int *pipvt = ipvt.fortran_vec ();
1251 retval = ComplexColumnVector (result, b_len); 1301
1252 } 1302 Array<Complex> z (nr);
1253 1303 Complex *pz = z.fortran_vec ();
1254 delete [] tmp_data; 1304
1255 delete [] ipvt; 1305 ComplexMatrix atmp = *this;
1256 delete [] z; 1306 Complex *tmp_data = atmp.fortran_vec ();
1307
1308 F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
1309
1310 if (f77_exception_encountered)
1311 (*current_liboctave_error_handler)
1312 ("unrecoverable error in dgeco");
1313 else
1314 {
1315 volatile double rcond_plus_one = rcond + 1.0;
1316
1317 if (rcond_plus_one == 1.0)
1318 {
1319 info = -2;
1320 }
1321 else
1322 {
1323 retval = b;
1324 Complex *result = retval.fortran_vec ();
1325
1326 F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt, result, 0));
1327
1328 if (f77_exception_encountered)
1329 (*current_liboctave_error_handler)
1330 ("unrecoverable error in dgesl");
1331 }
1332 }
1333 }
1257 1334
1258 return retval; 1335 return retval;
1259 } 1336 }
1260 1337
1261 ComplexMatrix 1338 ComplexMatrix
1274 } 1351 }
1275 1352
1276 ComplexMatrix 1353 ComplexMatrix
1277 ComplexMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const 1354 ComplexMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
1278 { 1355 {
1356 ComplexMatrix retval;
1357
1279 int nrhs = b.cols (); 1358 int nrhs = b.cols ();
1280 1359
1281 int m = rows (); 1360 int m = rows ();
1282 int n = cols (); 1361 int n = cols ();
1283 1362
1284 if (m == 0 || n == 0 || m != b.rows ()) 1363 if (m == 0 || n == 0 || m != b.rows ())
1285 { 1364 (*current_liboctave_error_handler)
1286 (*current_liboctave_error_handler) 1365 ("matrix dimension mismatch solution of linear equations");
1287 ("matrix dimension mismatch solution of linear equations");
1288 return Matrix ();
1289 }
1290
1291 Complex *tmp_data = dup (data (), length ());
1292
1293 int nrr = m > n ? m : n;
1294 ComplexMatrix result (nrr, nrhs);
1295
1296 for (int j = 0; j < nrhs; j++)
1297 for (int i = 0; i < m; i++)
1298 result.elem (i, j) = b.elem (i, j);
1299
1300 Complex *presult = result.fortran_vec ();
1301
1302 int len_s = m < n ? m : n;
1303 double *s = new double [len_s];
1304 double rcond = -1.0;
1305 int lwork;
1306 if (m < n)
1307 lwork = 2*m + (nrhs > n ? nrhs : n);
1308 else 1366 else
1309 lwork = 2*n + (nrhs > m ? nrhs : m); 1367 {
1310 1368 ComplexMatrix atmp = *this;
1311 Complex *work = new Complex [lwork]; 1369 Complex *tmp_data = atmp.fortran_vec ();
1312 1370
1313 int lrwork = (5 * (m < n ? m : n)) - 4; 1371 int nrr = m > n ? m : n;
1314 lrwork = lrwork > 1 ? lrwork : 1; 1372 ComplexMatrix result (nrr, nrhs);
1315 double *rwork = new double [lrwork]; 1373
1316 1374 for (int j = 0; j < nrhs; j++)
1317 F77_FCN (zgelss, ZGELSS) (m, n, nrhs, tmp_data, m, presult, nrr, s, 1375 for (int i = 0; i < m; i++)
1318 rcond, rank, work, lwork, rwork, info); 1376 result.elem (i, j) = b.elem (i, j);
1319 1377
1320 ComplexMatrix retval (n, nrhs); 1378 Complex *presult = result.fortran_vec ();
1321 for (int j = 0; j < nrhs; j++) 1379
1322 for (int i = 0; i < n; i++) 1380 int len_s = m < n ? m : n;
1323 retval.elem (i, j) = result.elem (i, j); 1381 Array<double> s (len_s);
1324 1382 double *ps = s.fortran_vec ();
1325 delete [] tmp_data; 1383 double rcond = -1.0;
1326 delete [] s; 1384 int lwork;
1327 delete [] work; 1385 if (m < n)
1328 delete [] rwork; 1386 lwork = 2*m + (nrhs > n ? nrhs : n);
1387 else
1388 lwork = 2*n + (nrhs > m ? nrhs : m);
1389
1390 Array<Complex> work (lwork);
1391 Complex *pwork = work.fortran_vec ();
1392
1393 int lrwork = (5 * (m < n ? m : n)) - 4;
1394 lrwork = lrwork > 1 ? lrwork : 1;
1395 Array<double> rwork (lrwork);
1396 double *prwork = rwork.fortran_vec ();
1397
1398 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
1399 nrr, ps, rcond, rank, pwork, lwork,
1400 prwork, info));
1401
1402 if (f77_exception_encountered)
1403 (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
1404 else
1405 {
1406 ComplexMatrix retval (n, nrhs);
1407 for (int j = 0; j < nrhs; j++)
1408 for (int i = 0; i < n; i++)
1409 retval.elem (i, j) = result.elem (i, j);
1410 }
1411 }
1329 1412
1330 return retval; 1413 return retval;
1331 } 1414 }
1332 1415
1333 ComplexColumnVector 1416 ComplexColumnVector
1347 1430
1348 ComplexColumnVector 1431 ComplexColumnVector
1349 ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info, 1432 ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info,
1350 int& rank) const 1433 int& rank) const
1351 { 1434 {
1435 ComplexColumnVector retval;
1436
1352 int nrhs = 1; 1437 int nrhs = 1;
1353 1438
1354 int m = rows (); 1439 int m = rows ();
1355 int n = cols (); 1440 int n = cols ();
1356 1441
1357 if (m == 0 || n == 0 || m != b.length ()) 1442 if (m == 0 || n == 0 || m != b.length ())
1358 { 1443 (*current_liboctave_error_handler)
1359 (*current_liboctave_error_handler) 1444 ("matrix dimension mismatch solution of least squares problem");
1360 ("matrix dimension mismatch solution of least squares problem");
1361 return ComplexColumnVector ();
1362 }
1363
1364 Complex *tmp_data = dup (data (), length ());
1365
1366 int nrr = m > n ? m : n;
1367 ComplexColumnVector result (nrr);
1368
1369 for (int i = 0; i < m; i++)
1370 result.elem (i) = b.elem (i);
1371
1372 Complex *presult = result.fortran_vec ();
1373
1374 int len_s = m < n ? m : n;
1375 double *s = new double [len_s];
1376 double rcond = -1.0;
1377 int lwork;
1378 if (m < n)
1379 lwork = 2*m + (nrhs > n ? nrhs : n);
1380 else 1445 else
1381 lwork = 2*n + (nrhs > m ? nrhs : m); 1446 {
1382 1447 ComplexMatrix atmp = *this;
1383 Complex *work = new Complex [lwork]; 1448 Complex *tmp_data = atmp.fortran_vec ();
1384 1449
1385 int lrwork = (5 * (m < n ? m : n)) - 4; 1450 int nrr = m > n ? m : n;
1386 lrwork = lrwork > 1 ? lrwork : 1; 1451 ComplexColumnVector result (nrr);
1387 double *rwork = new double [lrwork]; 1452
1388 1453 for (int i = 0; i < m; i++)
1389 F77_FCN (zgelss, ZGELSS) (m, n, nrhs, tmp_data, m, presult, nrr, s, 1454 result.elem (i) = b.elem (i);
1390 rcond, rank, work, lwork, rwork, info); 1455
1391 1456 Complex *presult = result.fortran_vec ();
1392 ComplexColumnVector retval (n); 1457
1393 for (int i = 0; i < n; i++) 1458 int len_s = m < n ? m : n;
1394 retval.elem (i) = result.elem (i); 1459 Array<double> s (len_s);
1395 1460 double *ps = s.fortran_vec ();
1396 delete [] tmp_data; 1461
1397 delete [] s; 1462 double rcond = -1.0;
1398 delete [] work; 1463
1399 delete [] rwork; 1464 int lwork;
1465 if (m < n)
1466 lwork = 2*m + (nrhs > n ? nrhs : n);
1467 else
1468 lwork = 2*n + (nrhs > m ? nrhs : m);
1469
1470 Array<Complex> work (lwork);
1471 Complex *pwork = work.fortran_vec ();
1472
1473 int lrwork = (5 * (m < n ? m : n)) - 4;
1474 lrwork = lrwork > 1 ? lrwork : 1;
1475 Array<double> rwork (lrwork);
1476 double *prwork = rwork.fortran_vec ();
1477
1478 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
1479 nrr, ps, rcond, rank, pwork, lwork,
1480 prwork, info));
1481
1482 if (f77_exception_encountered)
1483 (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
1484 else
1485 {
1486 ComplexColumnVector retval (n);
1487 for (int i = 0; i < n; i++)
1488 retval.elem (i) = result.elem (i);
1489 }
1490 }
1400 1491
1401 return retval; 1492 return retval;
1402 } 1493 }
1403 1494
1404 // Constants for matrix exponential calculation. 1495 // Constants for matrix exponential calculation.
1536 } 1627 }
1537 1628
1538 ComplexMatrix 1629 ComplexMatrix
1539 operator * (const ComplexColumnVector& v, const ComplexRowVector& a) 1630 operator * (const ComplexColumnVector& v, const ComplexRowVector& a)
1540 { 1631 {
1632 ComplexMatrix retval;
1633
1541 int len = v.length (); 1634 int len = v.length ();
1542 int a_len = a.length (); 1635 int a_len = a.length ();
1636
1543 if (len != a_len) 1637 if (len != a_len)
1544 { 1638 (*current_liboctave_error_handler)
1545 (*current_liboctave_error_handler) 1639 ("nonconformant vector multiplication attempted");
1546 ("nonconformant vector multiplication attempted"); 1640 else
1547 return ComplexMatrix (); 1641 {
1548 } 1642 if (len != 0)
1549 1643 {
1550 if (len == 0) 1644 retval.resize (len, a_len);
1551 return ComplexMatrix (len, len, 0.0); 1645 Complex *c = retval.fortran_vec ();
1552 1646
1553 Complex *c = new Complex [len * a_len]; 1647 F77_XFCN (zgemm, ZGEMM, ("N", "N", len, a_len, 1, 1.0,
1554 1648 v.data (), len, a.data (), 1, 0.0,
1555 F77_FCN (zgemm, ZGEMM) ("N", "N", len, a_len, 1, 1.0, v.data (), 1649 c, len, 1L, 1L));
1556 len, a.data (), 1, 0.0, c, len, 1L, 1L); 1650
1557 1651 if (f77_exception_encountered)
1558 return ComplexMatrix (c, len, a_len); 1652 (*current_liboctave_error_handler)
1653 ("unrecoverable error in zgemm");
1654 }
1655 }
1656
1657 return retval;
1559 } 1658 }
1560 1659
1561 // diagonal matrix by scalar -> matrix operations 1660 // diagonal matrix by scalar -> matrix operations
1562 1661
1563 ComplexMatrix 1662 ComplexMatrix
1764 return result; 1863 return result;
1765 } 1864 }
1766 1865
1767 ComplexMatrix 1866 ComplexMatrix
1768 operator * (const Matrix& m, const ComplexDiagMatrix& a) 1867 operator * (const Matrix& m, const ComplexDiagMatrix& a)
1868 {
1869 ComplexMatrix retval;
1870
1871 int nr = m.rows ();
1872 int nc = m.cols ();
1873
1874 int a_nr = a.rows ();
1875 int a_nc = a.cols ();
1876
1877 if (nc != a_nr)
1878 (*current_liboctave_error_handler)
1879 ("nonconformant matrix multiplication attempted");
1880 else
1881 {
1882 if (nr == 0 || nc == 0 || a_nc == 0)
1883 retval.resize (nr, a_nc, 0.0);
1884 else
1885 {
1886 retval.resize (nr, a_nc);
1887 Complex *c = retval.fortran_vec ();
1888
1889 Complex *ctmp = 0;
1890
1891 for (int j = 0; j < a.length (); j++)
1892 {
1893 int idx = j * nr;
1894 ctmp = c + idx;
1895 if (a.elem (j, j) == 1.0)
1896 {
1897 for (int i = 0; i < nr; i++)
1898 ctmp[i] = m.elem (i, j);
1899 }
1900 else if (a.elem (j, j) == 0.0)
1901 {
1902 for (int i = 0; i < nr; i++)
1903 ctmp[i] = 0.0;
1904 }
1905 else
1906 {
1907 for (int i = 0; i < nr; i++)
1908 ctmp[i] = a.elem (j, j) * m.elem (i, j);
1909 }
1910 }
1911
1912 if (a_nr < a_nc)
1913 {
1914 for (int i = nr * nc; i < nr * a_nc; i++)
1915 ctmp[i] = 0.0;
1916 }
1917 }
1918 }
1919
1920 return retval;
1921 }
1922
1923 // diagonal matrix by matrix -> matrix operations
1924
1925 ComplexMatrix
1926 operator + (const DiagMatrix& m, const ComplexMatrix& a)
1927 {
1928 int nr = m.rows ();
1929 int nc = m.cols ();
1930 if (nr != a.rows () || nc != a.cols ())
1931 {
1932 (*current_liboctave_error_handler)
1933 ("nonconformant matrix addition attempted");
1934 return ComplexMatrix ();
1935 }
1936
1937 if (nr == 0 || nc == 0)
1938 return ComplexMatrix (nr, nc);
1939
1940 ComplexMatrix result (a);
1941 for (int i = 0; i < m.length (); i++)
1942 result.elem (i, i) += m.elem (i, i);
1943
1944 return result;
1945 }
1946
1947 ComplexMatrix
1948 operator - (const DiagMatrix& m, const ComplexMatrix& a)
1949 {
1950 int nr = m.rows ();
1951 int nc = m.cols ();
1952 if (nr != a.rows () || nc != a.cols ())
1953 {
1954 (*current_liboctave_error_handler)
1955 ("nonconformant matrix subtraction attempted");
1956 return ComplexMatrix ();
1957 }
1958
1959 if (nr == 0 || nc == 0)
1960 return ComplexMatrix (nr, nc);
1961
1962 ComplexMatrix result (-a);
1963 for (int i = 0; i < m.length (); i++)
1964 result.elem (i, i) += m.elem (i, i);
1965
1966 return result;
1967 }
1968
1969 ComplexMatrix
1970 operator * (const DiagMatrix& m, const ComplexMatrix& a)
1769 { 1971 {
1770 int nr = m.rows (); 1972 int nr = m.rows ();
1771 int nc = m.cols (); 1973 int nc = m.cols ();
1772 int a_nr = a.rows (); 1974 int a_nr = a.rows ();
1773 int a_nc = a.cols (); 1975 int a_nc = a.cols ();
1777 ("nonconformant matrix multiplication attempted"); 1979 ("nonconformant matrix multiplication attempted");
1778 return ComplexMatrix (); 1980 return ComplexMatrix ();
1779 } 1981 }
1780 1982
1781 if (nr == 0 || nc == 0 || a_nc == 0) 1983 if (nr == 0 || nc == 0 || a_nc == 0)
1782 return ComplexMatrix (nr, a_nc, 0.0); 1984 return ComplexMatrix (nr, nc, 0.0);
1783 1985
1784 Complex *c = new Complex [nr*a_nc]; 1986 ComplexMatrix c (nr, a_nc);
1785 Complex *ctmp = 0; 1987
1786 1988 for (int i = 0; i < m.length (); i++)
1787 for (int j = 0; j < a.length (); j++) 1989 {
1788 { 1990 if (m.elem (i, i) == 1.0)
1789 int idx = j * nr; 1991 {
1790 ctmp = c + idx; 1992 for (int j = 0; j < a_nc; j++)
1791 if (a.elem (j, j) == 1.0) 1993 c.elem (i, j) = a.elem (i, j);
1792 { 1994 }
1793 for (int i = 0; i < nr; i++) 1995 else if (m.elem (i, i) == 0.0)
1794 ctmp[i] = m.elem (i, j); 1996 {
1795 } 1997 for (int j = 0; j < a_nc; j++)
1796 else if (a.elem (j, j) == 0.0) 1998 c.elem (i, j) = 0.0;
1797 {
1798 for (int i = 0; i < nr; i++)
1799 ctmp[i] = 0.0;
1800 } 1999 }
1801 else 2000 else
1802 { 2001 {
1803 for (int i = 0; i < nr; i++) 2002 for (int j = 0; j < a_nc; j++)
1804 ctmp[i] = a.elem (j, j) * m.elem (i, j); 2003 c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
1805 } 2004 }
1806 } 2005 }
1807 2006
1808 if (a_nr < a_nc) 2007 if (nr > nc)
1809 { 2008 {
1810 for (int i = nr * nc; i < nr * a_nc; i++) 2009 for (int j = 0; j < a_nc; j++)
1811 ctmp[i] = 0.0; 2010 for (int i = a_nr; i < nr; i++)
1812 } 2011 c.elem (i, j) = 0.0;
1813 2012 }
1814 return ComplexMatrix (c, nr, a_nc); 2013
1815 } 2014 return c;
1816 2015 }
1817 // diagonal matrix by matrix -> matrix operations 2016
1818 2017 ComplexMatrix
1819 ComplexMatrix 2018 operator + (const ComplexDiagMatrix& m, const Matrix& a)
1820 operator + (const DiagMatrix& m, const ComplexMatrix& a)
1821 { 2019 {
1822 int nr = m.rows (); 2020 int nr = m.rows ();
1823 int nc = m.cols (); 2021 int nc = m.cols ();
1824 if (nr != a.rows () || nc != a.cols ()) 2022 if (nr != a.rows () || nc != a.cols ())
1825 { 2023 {
1837 2035
1838 return result; 2036 return result;
1839 } 2037 }
1840 2038
1841 ComplexMatrix 2039 ComplexMatrix
1842 operator - (const DiagMatrix& m, const ComplexMatrix& a) 2040 operator - (const ComplexDiagMatrix& m, const Matrix& a)
1843 { 2041 {
1844 int nr = m.rows (); 2042 int nr = m.rows ();
1845 int nc = m.cols (); 2043 int nc = m.cols ();
1846 if (nr != a.rows () || nc != a.cols ()) 2044 if (nr != a.rows () || nc != a.cols ())
1847 { 2045 {
1859 2057
1860 return result; 2058 return result;
1861 } 2059 }
1862 2060
1863 ComplexMatrix 2061 ComplexMatrix
1864 operator * (const DiagMatrix& m, const ComplexMatrix& a) 2062 operator * (const ComplexDiagMatrix& m, const Matrix& a)
1865 { 2063 {
1866 int nr = m.rows (); 2064 int nr = m.rows ();
1867 int nc = m.cols (); 2065 int nc = m.cols ();
1868 int a_nr = a.rows (); 2066 int a_nr = a.rows ();
1869 int a_nc = a.cols (); 2067 int a_nc = a.cols ();
1873 ("nonconformant matrix multiplication attempted"); 2071 ("nonconformant matrix multiplication attempted");
1874 return ComplexMatrix (); 2072 return ComplexMatrix ();
1875 } 2073 }
1876 2074
1877 if (nr == 0 || nc == 0 || a_nc == 0) 2075 if (nr == 0 || nc == 0 || a_nc == 0)
1878 return ComplexMatrix (nr, nc, 0.0); 2076 return ComplexMatrix (nr, a_nc, 0.0);
1879 2077
1880 ComplexMatrix c (nr, a_nc); 2078 ComplexMatrix c (nr, a_nc);
1881 2079
1882 for (int i = 0; i < m.length (); i++) 2080 for (int i = 0; i < m.length (); i++)
1883 { 2081 {
1907 2105
1908 return c; 2106 return c;
1909 } 2107 }
1910 2108
1911 ComplexMatrix 2109 ComplexMatrix
1912 operator + (const ComplexDiagMatrix& m, const Matrix& a) 2110 operator + (const ComplexDiagMatrix& m, const ComplexMatrix& a)
1913 { 2111 {
1914 int nr = m.rows (); 2112 int nr = m.rows ();
1915 int nc = m.cols (); 2113 int nc = m.cols ();
1916 if (nr != a.rows () || nc != a.cols ()) 2114 if (nr != a.rows () || nc != a.cols ())
1917 { 2115 {
1929 2127
1930 return result; 2128 return result;
1931 } 2129 }
1932 2130
1933 ComplexMatrix 2131 ComplexMatrix
1934 operator - (const ComplexDiagMatrix& m, const Matrix& a) 2132 operator - (const ComplexDiagMatrix& m, const ComplexMatrix& a)
1935 { 2133 {
1936 int nr = m.rows (); 2134 int nr = m.rows ();
1937 int nc = m.cols (); 2135 int nc = m.cols ();
1938 if (nr != a.rows () || nc != a.cols ()) 2136 if (nr != a.rows () || nc != a.cols ())
1939 { 2137 {
1951 2149
1952 return result; 2150 return result;
1953 } 2151 }
1954 2152
1955 ComplexMatrix 2153 ComplexMatrix
1956 operator * (const ComplexDiagMatrix& m, const Matrix& a) 2154 operator * (const ComplexDiagMatrix& m, const ComplexMatrix& a)
1957 { 2155 {
1958 int nr = m.rows (); 2156 int nr = m.rows ();
1959 int nc = m.cols (); 2157 int nc = m.cols ();
1960 int a_nr = a.rows (); 2158 int a_nr = a.rows ();
1961 int a_nc = a.cols (); 2159 int a_nc = a.cols ();
1998 } 2196 }
1999 2197
2000 return c; 2198 return c;
2001 } 2199 }
2002 2200
2003 ComplexMatrix 2201 // matrix by matrix -> matrix operations
2004 operator + (const ComplexDiagMatrix& m, const ComplexMatrix& a) 2202
2203 ComplexMatrix&
2204 ComplexMatrix::operator += (const Matrix& a)
2205 {
2206 int nr = rows ();
2207 int nc = cols ();
2208 if (nr != a.rows () || nc != a.cols ())
2209 {
2210 (*current_liboctave_error_handler)
2211 ("nonconformant matrix += operation attempted");
2212 return *this;
2213 }
2214
2215 if (nr == 0 || nc == 0)
2216 return *this;
2217
2218 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
2219
2220 add2 (d, a.data (), length ());
2221 return *this;
2222 }
2223
2224 ComplexMatrix&
2225 ComplexMatrix::operator -= (const Matrix& a)
2226 {
2227 int nr = rows ();
2228 int nc = cols ();
2229 if (nr != a.rows () || nc != a.cols ())
2230 {
2231 (*current_liboctave_error_handler)
2232 ("nonconformant matrix -= operation attempted");
2233 return *this;
2234 }
2235
2236 if (nr == 0 || nc == 0)
2237 return *this;
2238
2239 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
2240
2241 subtract2 (d, a.data (), length ());
2242 return *this;
2243 }
2244
2245 ComplexMatrix&
2246 ComplexMatrix::operator += (const ComplexMatrix& a)
2247 {
2248 int nr = rows ();
2249 int nc = cols ();
2250 if (nr != a.rows () || nc != a.cols ())
2251 {
2252 (*current_liboctave_error_handler)
2253 ("nonconformant matrix += operation attempted");
2254 return *this;
2255 }
2256
2257 if (nr == 0 || nc == 0)
2258 return *this;
2259
2260 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
2261
2262 add2 (d, a.data (), length ());
2263 return *this;
2264 }
2265
2266 ComplexMatrix&
2267 ComplexMatrix::operator -= (const ComplexMatrix& a)
2268 {
2269 int nr = rows ();
2270 int nc = cols ();
2271 if (nr != a.rows () || nc != a.cols ())
2272 {
2273 (*current_liboctave_error_handler)
2274 ("nonconformant matrix -= operation attempted");
2275 return *this;
2276 }
2277
2278 if (nr == 0 || nc == 0)
2279 return *this;
2280
2281 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
2282
2283 subtract2 (d, a.data (), length ());
2284 return *this;
2285 }
2286
2287 // unary operations
2288
2289 Matrix
2290 ComplexMatrix::operator ! (void) const
2291 {
2292 return Matrix (not (data (), length ()), rows (), cols ());
2293 }
2294
2295 // matrix by scalar -> matrix operations
2296
2297 ComplexMatrix
2298 operator + (const Matrix& a, const Complex& s)
2299 {
2300 return ComplexMatrix (add (a.data (), a.length (), s),
2301 a.rows (), a.cols ());
2302 }
2303
2304 ComplexMatrix
2305 operator - (const Matrix& a, const Complex& s)
2306 {
2307 return ComplexMatrix (subtract (a.data (), a.length (), s),
2308 a.rows (), a.cols ());
2309 }
2310
2311 ComplexMatrix
2312 operator * (const Matrix& a, const Complex& s)
2313 {
2314 return ComplexMatrix (multiply (a.data (), a.length (), s),
2315 a.rows (), a.cols ());
2316 }
2317
2318 ComplexMatrix
2319 operator / (const Matrix& a, const Complex& s)
2320 {
2321 return ComplexMatrix (divide (a.data (), a.length (), s),
2322 a.rows (), a.cols ());
2323 }
2324
2325 ComplexMatrix
2326 operator + (const ComplexMatrix& a, double s)
2327 {
2328 return ComplexMatrix (add (a.data (), a.length (), s),
2329 a.rows (), a.cols ());
2330 }
2331
2332 ComplexMatrix
2333 operator - (const ComplexMatrix& a, double s)
2334 {
2335 return ComplexMatrix (subtract (a.data (), a.length (), s),
2336 a.rows (), a.cols ());
2337 }
2338
2339 ComplexMatrix
2340 operator * (const ComplexMatrix& a, double s)
2341 {
2342 return ComplexMatrix (multiply (a.data (), a.length (), s),
2343 a.rows (), a.cols ());
2344 }
2345
2346 ComplexMatrix
2347 operator / (const ComplexMatrix& a, double s)
2348 {
2349 return ComplexMatrix (divide (a.data (), a.length (), s),
2350 a.rows (), a.cols ());
2351 }
2352
2353 // scalar by matrix -> matrix operations
2354
2355 ComplexMatrix
2356 operator + (double s, const ComplexMatrix& a)
2357 {
2358 return ComplexMatrix (add (a.data (), a.length (), s), a.rows (),
2359 a.cols ());
2360 }
2361
2362 ComplexMatrix
2363 operator - (double s, const ComplexMatrix& a)
2364 {
2365 return ComplexMatrix (subtract (s, a.data (), a.length ()),
2366 a.rows (), a.cols ());
2367 }
2368
2369 ComplexMatrix
2370 operator * (double s, const ComplexMatrix& a)
2371 {
2372 return ComplexMatrix (multiply (a.data (), a.length (), s),
2373 a.rows (), a.cols ());
2374 }
2375
2376 ComplexMatrix
2377 operator / (double s, const ComplexMatrix& a)
2378 {
2379 return ComplexMatrix (divide (s, a.data (), a.length ()),
2380 a.rows (), a.cols ());
2381 }
2382
2383 ComplexMatrix
2384 operator + (const Complex& s, const Matrix& a)
2385 {
2386 return ComplexMatrix (add (s, a.data (), a.length ()),
2387 a.rows (), a.cols ());
2388 }
2389
2390 ComplexMatrix
2391 operator - (const Complex& s, const Matrix& a)
2392 {
2393 return ComplexMatrix (subtract (s, a.data (), a.length ()),
2394 a.rows (), a.cols ());
2395 }
2396
2397 ComplexMatrix
2398 operator * (const Complex& s, const Matrix& a)
2399 {
2400 return ComplexMatrix (multiply (a.data (), a.length (), s),
2401 a.rows (), a.cols ());
2402 }
2403
2404 ComplexMatrix
2405 operator / (const Complex& s, const Matrix& a)
2406 {
2407 return ComplexMatrix (divide (s, a.data (), a.length ()),
2408 a.rows (), a.cols ());
2409 }
2410
2411 // matrix by diagonal matrix -> matrix operations
2412
2413 ComplexMatrix
2414 operator + (const ComplexMatrix& m, const DiagMatrix& a)
2005 { 2415 {
2006 int nr = m.rows (); 2416 int nr = m.rows ();
2007 int nc = m.cols (); 2417 int nc = m.cols ();
2008 if (nr != a.rows () || nc != a.cols ()) 2418 if (nr != a.rows () || nc != a.cols ())
2009 { 2419 {
2013 } 2423 }
2014 2424
2015 if (nr == 0 || nc == 0) 2425 if (nr == 0 || nc == 0)
2016 return ComplexMatrix (nr, nc); 2426 return ComplexMatrix (nr, nc);
2017 2427
2018 ComplexMatrix result (a); 2428 ComplexMatrix result (m);
2019 for (int i = 0; i < m.length (); i++) 2429 for (int i = 0; i < a.length (); i++)
2020 result.elem (i, i) += m.elem (i, i); 2430 result.elem (i, i) += a.elem (i, i);
2021 2431
2022 return result; 2432 return result;
2023 } 2433 }
2024 2434
2025 ComplexMatrix 2435 ComplexMatrix
2026 operator - (const ComplexDiagMatrix& m, const ComplexMatrix& a) 2436 operator - (const ComplexMatrix& m, const DiagMatrix& a)
2027 { 2437 {
2028 int nr = m.rows (); 2438 int nr = m.rows ();
2029 int nc = m.cols (); 2439 int nc = m.cols ();
2030 if (nr != a.rows () || nc != a.cols ()) 2440 if (nr != a.rows () || nc != a.cols ())
2031 { 2441 {
2035 } 2445 }
2036 2446
2037 if (nr == 0 || nc == 0) 2447 if (nr == 0 || nc == 0)
2038 return ComplexMatrix (nr, nc); 2448 return ComplexMatrix (nr, nc);
2039 2449
2040 ComplexMatrix result (-a); 2450 ComplexMatrix result (m);
2041 for (int i = 0; i < m.length (); i++) 2451 for (int i = 0; i < a.length (); i++)
2042 result.elem (i, i) += m.elem (i, i); 2452 result.elem (i, i) -= a.elem (i, i);
2043 2453
2044 return result; 2454 return result;
2045 } 2455 }
2046 2456
2047 ComplexMatrix 2457 ComplexMatrix
2048 operator * (const ComplexDiagMatrix& m, const ComplexMatrix& a) 2458 operator * (const ComplexMatrix& m, const DiagMatrix& a)
2049 { 2459 {
2460 ComplexMatrix retval;
2461
2050 int nr = m.rows (); 2462 int nr = m.rows ();
2051 int nc = m.cols (); 2463 int nc = m.cols ();
2052 int a_nr = a.rows (); 2464
2053 int a_nc = a.cols (); 2465 int a_nc = a.cols ();
2054 if (nc != a_nr) 2466
2055 { 2467 if (nc != a.rows ())
2056 (*current_liboctave_error_handler) 2468 (*current_liboctave_error_handler)
2057 ("nonconformant matrix multiplication attempted"); 2469 ("nonconformant matrix multiplication attempted");
2058 return ComplexMatrix (); 2470 else
2059 } 2471 {
2060 2472 if (nr == 0 || nc == 0 || a_nc == 0)
2061 if (nr == 0 || nc == 0 || a_nc == 0) 2473 retval.resize (nr, nc, 0.0);
2062 return ComplexMatrix (nr, a_nc, 0.0);
2063
2064 ComplexMatrix c (nr, a_nc);
2065
2066 for (int i = 0; i < m.length (); i++)
2067 {
2068 if (m.elem (i, i) == 1.0)
2069 {
2070 for (int j = 0; j < a_nc; j++)
2071 c.elem (i, j) = a.elem (i, j);
2072 }
2073 else if (m.elem (i, i) == 0.0)
2074 {
2075 for (int j = 0; j < a_nc; j++)
2076 c.elem (i, j) = 0.0;
2077 }
2078 else 2474 else
2079 { 2475 {
2080 for (int j = 0; j < a_nc; j++) 2476 retval.resize (nr, a_nc);
2081 c.elem (i, j) = m.elem (i, i) * a.elem (i, j); 2477 Complex *c = retval.fortran_vec ();
2082 } 2478 Complex *ctmp = 0;
2083 } 2479
2084 2480 for (int j = 0; j < a.length (); j++)
2085 if (nr > nc) 2481 {
2086 { 2482 int idx = j * nr;
2087 for (int j = 0; j < a_nc; j++) 2483 ctmp = c + idx;
2088 for (int i = a_nr; i < nr; i++) 2484 if (a.elem (j, j) == 1.0)
2089 c.elem (i, j) = 0.0; 2485 {
2090 } 2486 for (int i = 0; i < nr; i++)
2091 2487 ctmp[i] = m.elem (i, j);
2092 return c; 2488 }
2093 } 2489 else if (a.elem (j, j) == 0.0)
2094 2490 {
2095 // matrix by matrix -> matrix operations 2491 for (int i = 0; i < nr; i++)
2096 2492 ctmp[i] = 0.0;
2097 ComplexMatrix& 2493 }
2098 ComplexMatrix::operator += (const Matrix& a) 2494 else
2099 { 2495 {
2100 int nr = rows (); 2496 for (int i = 0; i < nr; i++)
2101 int nc = cols (); 2497 ctmp[i] = a.elem (j, j) * m.elem (i, j);
2102 if (nr != a.rows () || nc != a.cols ()) 2498 }
2103 { 2499 }
2104 (*current_liboctave_error_handler) 2500
2105 ("nonconformant matrix += operation attempted"); 2501 if (a.rows () < a_nc)
2106 return *this; 2502 {
2107 } 2503 for (int i = nr * nc; i < nr * a_nc; i++)
2108 2504 ctmp[i] = 0.0;
2109 if (nr == 0 || nc == 0) 2505 }
2110 return *this; 2506 }
2111 2507 }
2112 Complex *d = fortran_vec (); // Ensures only one reference to my privates! 2508
2113 2509 return retval;
2114 add2 (d, a.data (), length ()); 2510 }
2115 return *this; 2511
2116 } 2512 ComplexMatrix
2117 2513 operator + (const ComplexMatrix& m, const ComplexDiagMatrix& a)
2118 ComplexMatrix&
2119 ComplexMatrix::operator -= (const Matrix& a)
2120 {
2121 int nr = rows ();
2122 int nc = cols ();
2123 if (nr != a.rows () || nc != a.cols ())
2124 {
2125 (*current_liboctave_error_handler)
2126 ("nonconformant matrix -= operation attempted");
2127 return *this;
2128 }
2129
2130 if (nr == 0 || nc == 0)
2131 return *this;
2132
2133 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
2134
2135 subtract2 (d, a.data (), length ());
2136 return *this;
2137 }
2138
2139 ComplexMatrix&
2140 ComplexMatrix::operator += (const ComplexMatrix& a)
2141 {
2142 int nr = rows ();
2143 int nc = cols ();
2144 if (nr != a.rows () || nc != a.cols ())
2145 {
2146 (*current_liboctave_error_handler)
2147 ("nonconformant matrix += operation attempted");
2148 return *this;
2149 }
2150
2151 if (nr == 0 || nc == 0)
2152 return *this;
2153
2154 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
2155
2156 add2 (d, a.data (), length ());
2157 return *this;
2158 }
2159
2160 ComplexMatrix&
2161 ComplexMatrix::operator -= (const ComplexMatrix& a)
2162 {
2163 int nr = rows ();
2164 int nc = cols ();
2165 if (nr != a.rows () || nc != a.cols ())
2166 {
2167 (*current_liboctave_error_handler)
2168 ("nonconformant matrix -= operation attempted");
2169 return *this;
2170 }
2171
2172 if (nr == 0 || nc == 0)
2173 return *this;
2174
2175 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
2176
2177 subtract2 (d, a.data (), length ());
2178 return *this;
2179 }
2180
2181 // unary operations
2182
2183 Matrix
2184 ComplexMatrix::operator ! (void) const
2185 {
2186 return Matrix (not (data (), length ()), rows (), cols ());
2187 }
2188
2189 // matrix by scalar -> matrix operations
2190
2191 ComplexMatrix
2192 operator + (const Matrix& a, const Complex& s)
2193 {
2194 return ComplexMatrix (add (a.data (), a.length (), s),
2195 a.rows (), a.cols ());
2196 }
2197
2198 ComplexMatrix
2199 operator - (const Matrix& a, const Complex& s)
2200 {
2201 return ComplexMatrix (subtract (a.data (), a.length (), s),
2202 a.rows (), a.cols ());
2203 }
2204
2205 ComplexMatrix
2206 operator * (const Matrix& a, const Complex& s)
2207 {
2208 return ComplexMatrix (multiply (a.data (), a.length (), s),
2209 a.rows (), a.cols ());
2210 }
2211
2212 ComplexMatrix
2213 operator / (const Matrix& a, const Complex& s)
2214 {
2215 return ComplexMatrix (divide (a.data (), a.length (), s),
2216 a.rows (), a.cols ());
2217 }
2218
2219 ComplexMatrix
2220 operator + (const ComplexMatrix& a, double s)
2221 {
2222 return ComplexMatrix (add (a.data (), a.length (), s),
2223 a.rows (), a.cols ());
2224 }
2225
2226 ComplexMatrix
2227 operator - (const ComplexMatrix& a, double s)
2228 {
2229 return ComplexMatrix (subtract (a.data (), a.length (), s),
2230 a.rows (), a.cols ());
2231 }
2232
2233 ComplexMatrix
2234 operator * (const ComplexMatrix& a, double s)
2235 {
2236 return ComplexMatrix (multiply (a.data (), a.length (), s),
2237 a.rows (), a.cols ());
2238 }
2239
2240 ComplexMatrix
2241 operator / (const ComplexMatrix& a, double s)
2242 {
2243 return ComplexMatrix (divide (a.data (), a.length (), s),
2244 a.rows (), a.cols ());
2245 }
2246
2247 // scalar by matrix -> matrix operations
2248
2249 ComplexMatrix
2250 operator + (double s, const ComplexMatrix& a)
2251 {
2252 return ComplexMatrix (add (a.data (), a.length (), s), a.rows (),
2253 a.cols ());
2254 }
2255
2256 ComplexMatrix
2257 operator - (double s, const ComplexMatrix& a)
2258 {
2259 return ComplexMatrix (subtract (s, a.data (), a.length ()),
2260 a.rows (), a.cols ());
2261 }
2262
2263 ComplexMatrix
2264 operator * (double s, const ComplexMatrix& a)
2265 {
2266 return ComplexMatrix (multiply (a.data (), a.length (), s),
2267 a.rows (), a.cols ());
2268 }
2269
2270 ComplexMatrix
2271 operator / (double s, const ComplexMatrix& a)
2272 {
2273 return ComplexMatrix (divide (s, a.data (), a.length ()),
2274 a.rows (), a.cols ());
2275 }
2276
2277 ComplexMatrix
2278 operator + (const Complex& s, const Matrix& a)
2279 {
2280 return ComplexMatrix (add (s, a.data (), a.length ()),
2281 a.rows (), a.cols ());
2282 }
2283
2284 ComplexMatrix
2285 operator - (const Complex& s, const Matrix& a)
2286 {
2287 return ComplexMatrix (subtract (s, a.data (), a.length ()),
2288 a.rows (), a.cols ());
2289 }
2290
2291 ComplexMatrix
2292 operator * (const Complex& s, const Matrix& a)
2293 {
2294 return ComplexMatrix (multiply (a.data (), a.length (), s),
2295 a.rows (), a.cols ());
2296 }
2297
2298 ComplexMatrix
2299 operator / (const Complex& s, const Matrix& a)
2300 {
2301 return ComplexMatrix (divide (s, a.data (), a.length ()),
2302 a.rows (), a.cols ());
2303 }
2304
2305 // matrix by diagonal matrix -> matrix operations
2306
2307 ComplexMatrix
2308 operator + (const ComplexMatrix& m, const DiagMatrix& a)
2309 { 2514 {
2310 int nr = m.rows (); 2515 int nr = m.rows ();
2311 int nc = m.cols (); 2516 int nc = m.cols ();
2312 if (nr != a.rows () || nc != a.cols ()) 2517 if (nr != a.rows () || nc != a.cols ())
2313 { 2518 {
2325 2530
2326 return result; 2531 return result;
2327 } 2532 }
2328 2533
2329 ComplexMatrix 2534 ComplexMatrix
2330 operator - (const ComplexMatrix& m, const DiagMatrix& a) 2535 operator - (const ComplexMatrix& m, const ComplexDiagMatrix& a)
2331 { 2536 {
2332 int nr = m.rows (); 2537 int nr = m.rows ();
2333 int nc = m.cols (); 2538 int nc = m.cols ();
2334 if (nr != a.rows () || nc != a.cols ()) 2539 if (nr != a.rows () || nc != a.cols ())
2335 { 2540 {
2347 2552
2348 return result; 2553 return result;
2349 } 2554 }
2350 2555
2351 ComplexMatrix 2556 ComplexMatrix
2352 operator * (const ComplexMatrix& m, const DiagMatrix& a) 2557 operator * (const ComplexMatrix& m, const ComplexDiagMatrix& a)
2353 { 2558 {
2559 ComplexMatrix retval;
2560
2354 int nr = m.rows (); 2561 int nr = m.rows ();
2355 int nc = m.cols (); 2562 int nc = m.cols ();
2563
2356 int a_nc = a.cols (); 2564 int a_nc = a.cols ();
2565
2357 if (nc != a.rows ()) 2566 if (nc != a.rows ())
2358 { 2567 (*current_liboctave_error_handler)
2359 (*current_liboctave_error_handler) 2568 ("nonconformant matrix multiplication attempted");
2360 ("nonconformant matrix multiplication attempted"); 2569 else
2361 return ComplexMatrix (); 2570 {
2362 } 2571 if (nr == 0 || nc == 0 || a_nc == 0)
2363 2572 retval.resize (nr, nc, 0.0);
2364 if (nr == 0 || nc == 0 || a_nc == 0)
2365 return ComplexMatrix (nr, nc, 0.0);
2366
2367 Complex *c = new Complex [nr*a_nc];
2368 Complex *ctmp = 0;
2369
2370 for (int j = 0; j < a.length (); j++)
2371 {
2372 int idx = j * nr;
2373 ctmp = c + idx;
2374 if (a.elem (j, j) == 1.0)
2375 {
2376 for (int i = 0; i < nr; i++)
2377 ctmp[i] = m.elem (i, j);
2378 }
2379 else if (a.elem (j, j) == 0.0)
2380 {
2381 for (int i = 0; i < nr; i++)
2382 ctmp[i] = 0.0;
2383 }
2384 else 2573 else
2385 { 2574 {
2386 for (int i = 0; i < nr; i++) 2575 retval.resize (nr, nc);
2387 ctmp[i] = a.elem (j, j) * m.elem (i, j); 2576 Complex *c = retval.fortran_vec ();
2388 } 2577 Complex *ctmp = 0;
2389 } 2578
2390 2579 for (int j = 0; j < a.length (); j++)
2391 if (a.rows () < a_nc) 2580 {
2392 { 2581 int idx = j * nr;
2393 for (int i = nr * nc; i < nr * a_nc; i++) 2582 ctmp = c + idx;
2394 ctmp[i] = 0.0; 2583 if (a.elem (j, j) == 1.0)
2395 } 2584 {
2396 2585 for (int i = 0; i < nr; i++)
2397 return ComplexMatrix (c, nr, a_nc); 2586 ctmp[i] = m.elem (i, j);
2398 } 2587 }
2399 2588 else if (a.elem (j, j) == 0.0)
2400 ComplexMatrix 2589 {
2401 operator + (const ComplexMatrix& m, const ComplexDiagMatrix& a) 2590 for (int i = 0; i < nr; i++)
2591 ctmp[i] = 0.0;
2592 }
2593 else
2594 {
2595 for (int i = 0; i < nr; i++)
2596 ctmp[i] = a.elem (j, j) * m.elem (i, j);
2597 }
2598 }
2599
2600 if (a.rows () < a_nc)
2601 {
2602 for (int i = nr * nc; i < nr * a_nc; i++)
2603 ctmp[i] = 0.0;
2604 }
2605 }
2606 }
2607
2608 return retval;
2609 }
2610
2611 // matrix by matrix -> matrix operations
2612
2613 ComplexMatrix
2614 operator + (const ComplexMatrix& m, const Matrix& a)
2402 { 2615 {
2403 int nr = m.rows (); 2616 int nr = m.rows ();
2404 int nc = m.cols (); 2617 int nc = m.cols ();
2405 if (nr != a.rows () || nc != a.cols ()) 2618 if (nr != a.rows () || nc != a.cols ())
2406 { 2619 {
2410 } 2623 }
2411 2624
2412 if (nr == 0 || nc == 0) 2625 if (nr == 0 || nc == 0)
2413 return ComplexMatrix (nr, nc); 2626 return ComplexMatrix (nr, nc);
2414 2627
2415 ComplexMatrix result (m); 2628 return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
2416 for (int i = 0; i < a.length (); i++) 2629 }
2417 result.elem (i, i) += a.elem (i, i); 2630
2418 2631 ComplexMatrix
2419 return result; 2632 operator - (const ComplexMatrix& m, const Matrix& a)
2420 }
2421
2422 ComplexMatrix
2423 operator - (const ComplexMatrix& m, const ComplexDiagMatrix& a)
2424 { 2633 {
2425 int nr = m.rows (); 2634 int nr = m.rows ();
2426 int nc = m.cols (); 2635 int nc = m.cols ();
2427 if (nr != a.rows () || nc != a.cols ()) 2636 if (nr != a.rows () || nc != a.cols ())
2428 { 2637 {
2432 } 2641 }
2433 2642
2434 if (nr == 0 || nc == 0) 2643 if (nr == 0 || nc == 0)
2435 return ComplexMatrix (nr, nc); 2644 return ComplexMatrix (nr, nc);
2436 2645
2437 ComplexMatrix result (m); 2646 return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc);
2438 for (int i = 0; i < a.length (); i++) 2647 }
2439 result.elem (i, i) -= a.elem (i, i); 2648
2440 2649 ComplexMatrix
2441 return result; 2650 operator + (const Matrix& m, const ComplexMatrix& a)
2442 }
2443
2444 ComplexMatrix
2445 operator * (const ComplexMatrix& m, const ComplexDiagMatrix& a)
2446 {
2447 int nr = m.rows ();
2448 int nc = m.cols ();
2449 int a_nc = a.cols ();
2450 if (nc != a.rows ())
2451 {
2452 (*current_liboctave_error_handler)
2453 ("nonconformant matrix multiplication attempted");
2454 return ComplexMatrix ();
2455 }
2456
2457 if (nr == 0 || nc == 0 || a_nc == 0)
2458 return ComplexMatrix (nr, nc, 0.0);
2459
2460 Complex *c = new Complex [nr*a_nc];
2461 Complex *ctmp = 0;
2462
2463 for (int j = 0; j < a.length (); j++)
2464 {
2465 int idx = j * nr;
2466 ctmp = c + idx;
2467 if (a.elem (j, j) == 1.0)
2468 {
2469 for (int i = 0; i < nr; i++)
2470 ctmp[i] = m.elem (i, j);
2471 }
2472 else if (a.elem (j, j) == 0.0)
2473 {
2474 for (int i = 0; i < nr; i++)
2475 ctmp[i] = 0.0;
2476 }
2477 else
2478 {
2479 for (int i = 0; i < nr; i++)
2480 ctmp[i] = a.elem (j, j) * m.elem (i, j);
2481 }
2482 }
2483
2484 if (a.rows () < a_nc)
2485 {
2486 for (int i = nr * nc; i < nr * a_nc; i++)
2487 ctmp[i] = 0.0;
2488 }
2489
2490 return ComplexMatrix (c, nr, a_nc);
2491 }
2492
2493 // matrix by matrix -> matrix operations
2494
2495 ComplexMatrix
2496 operator + (const ComplexMatrix& m, const Matrix& a)
2497 { 2651 {
2498 int nr = m.rows (); 2652 int nr = m.rows ();
2499 int nc = m.cols (); 2653 int nc = m.cols ();
2500 if (nr != a.rows () || nc != a.cols ()) 2654 if (nr != a.rows () || nc != a.cols ())
2501 { 2655 {
2502 (*current_liboctave_error_handler) 2656 (*current_liboctave_error_handler)
2503 ("nonconformant matrix addition attempted"); 2657 ("nonconformant matrix addition attempted");
2504 return ComplexMatrix (); 2658 return ComplexMatrix ();
2505 } 2659 }
2506 2660
2507 if (nr == 0 || nc == 0)
2508 return ComplexMatrix (nr, nc);
2509
2510 return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc); 2661 return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
2511 } 2662 }
2512 2663
2513 ComplexMatrix 2664 ComplexMatrix
2514 operator - (const ComplexMatrix& m, const Matrix& a) 2665 operator - (const Matrix& m, const ComplexMatrix& a)
2515 { 2666 {
2516 int nr = m.rows (); 2667 int nr = m.rows ();
2517 int nc = m.cols (); 2668 int nc = m.cols ();
2518 if (nr != a.rows () || nc != a.cols ()) 2669 if (nr != a.rows () || nc != a.cols ())
2519 { 2670 {
2527 2678
2528 return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc); 2679 return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc);
2529 } 2680 }
2530 2681
2531 ComplexMatrix 2682 ComplexMatrix
2532 operator + (const Matrix& m, const ComplexMatrix& a) 2683 operator * (const ComplexMatrix& m, const Matrix& a)
2533 { 2684 {
2685 ComplexMatrix tmp (a);
2686 return m * tmp;
2687 }
2688
2689 ComplexMatrix
2690 operator * (const Matrix& m, const ComplexMatrix& a)
2691 {
2692 ComplexMatrix tmp (m);
2693 return tmp * a;
2694 }
2695
2696 ComplexMatrix
2697 operator * (const ComplexMatrix& m, const ComplexMatrix& a)
2698 {
2699 ComplexMatrix retval;
2700
2534 int nr = m.rows (); 2701 int nr = m.rows ();
2535 int nc = m.cols (); 2702 int nc = m.cols ();
2536 if (nr != a.rows () || nc != a.cols ()) 2703
2537 {
2538 (*current_liboctave_error_handler)
2539 ("nonconformant matrix addition attempted");
2540 return ComplexMatrix ();
2541 }
2542
2543 return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
2544 }
2545
2546 ComplexMatrix
2547 operator - (const Matrix& m, const ComplexMatrix& a)
2548 {
2549 int nr = m.rows ();
2550 int nc = m.cols ();
2551 if (nr != a.rows () || nc != a.cols ())
2552 {
2553 (*current_liboctave_error_handler)
2554 ("nonconformant matrix subtraction attempted");
2555 return ComplexMatrix ();
2556 }
2557
2558 if (nr == 0 || nc == 0)
2559 return ComplexMatrix (nr, nc);
2560
2561 return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc);
2562 }
2563
2564 ComplexMatrix
2565 operator * (const ComplexMatrix& m, const Matrix& a)
2566 {
2567 ComplexMatrix tmp (a);
2568 return m * tmp;
2569 }
2570
2571 ComplexMatrix
2572 operator * (const Matrix& m, const ComplexMatrix& a)
2573 {
2574 ComplexMatrix tmp (m);
2575 return tmp * a;
2576 }
2577
2578 ComplexMatrix
2579 operator * (const ComplexMatrix& m, const ComplexMatrix& a)
2580 {
2581 int nr = m.rows ();
2582 int nc = m.cols ();
2583 int a_nc = a.cols (); 2704 int a_nc = a.cols ();
2705
2584 if (nc != a.rows ()) 2706 if (nc != a.rows ())
2585 { 2707 (*current_liboctave_error_handler)
2586 (*current_liboctave_error_handler) 2708 ("nonconformant matrix multiplication attempted");
2587 ("nonconformant matrix multiplication attempted"); 2709 else
2588 return ComplexMatrix (); 2710 {
2589 } 2711 if (nr == 0 || nc == 0 || a_nc == 0)
2590 2712 retval.resize (nr, nc, 0.0);
2591 if (nr == 0 || nc == 0 || a_nc == 0) 2713 else
2592 return ComplexMatrix (nr, nc, 0.0); 2714 {
2593 2715 int ld = nr;
2594 int ld = nr; 2716 int lda = a.rows ();
2595 int lda = a.rows (); 2717
2596 2718 retval.resize (nr, a_nc);
2597 Complex *c = new Complex [nr*a_nc]; 2719 Complex *c = retval.fortran_vec ();
2598 2720
2599 F77_FCN (zgemm, ZGEMM) ("N", "N", nr, a_nc, nc, 1.0, m.data (), ld, 2721 F77_XFCN (zgemm, ZGEMM, ("N", "N", nr, a_nc, nc, 1.0,
2600 a.data (), lda, 0.0, c, nr, 1L, 1L); 2722 m.data (), ld, a.data (), lda, 0.0,
2601 2723 c, nr, 1L, 1L));
2602 return ComplexMatrix (c, nr, a_nc); 2724
2725 if (f77_exception_encountered)
2726 (*current_liboctave_error_handler)
2727 ("unrecoverable error in zgemm");
2728 }
2729 }
2730
2731 return retval;
2603 } 2732 }
2604 2733
2605 ComplexMatrix 2734 ComplexMatrix
2606 product (const ComplexMatrix& m, const Matrix& a) 2735 product (const ComplexMatrix& m, const Matrix& a)
2607 { 2736 {