/[escript]/trunk/escript/src/DataArrayView.h
ViewVC logotype

Diff of /trunk/escript/src/DataArrayView.h

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 757 by woo409, Mon Jun 26 13:12:56 2006 UTC revision 775 by ksteube, Mon Jul 10 04:00:08 2006 UTC
# Line 18  Line 18 
18    
19  #include "esysUtils/EsysAssert.h"  #include "esysUtils/EsysAssert.h"
20    
21    #include "DataException.h"
22  #include "DataVector.h"  #include "DataVector.h"
23  #include "LocalOps.h"  #include "LocalOps.h"
24    
# Line 861  class DataArrayView { Line 862  class DataArrayView {
862    
863    /**    /**
864       \brief       \brief
865         computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2
866    
867         \param in - Input - matrix
868         \param inOffset - Input - offset into in
869         \param ev - Output - The symmetric matrix
870         \param inOffset - Input - offset into ev
871      */
872      static
873      inline
874      void
875      symmetric(DataArrayView& in,
876                ValueType::size_type inOffset,
877                DataArrayView& ev,
878                ValueType::size_type evOffset)
879      {
880       if (in.getRank() == 2) {
881         int i0, i1;
882         int s0=in.getShape()[0];
883         int s1=in.getShape()[1];
884         for (i0=0; i0<s0; i0++) {
885           for (i1=0; i1<s1; i1++) {
886             (*(ev.m_data))[evOffset+ev.index(i0,i1)] = ((*(in.m_data))[inOffset+in.index(i0,i1)] + (*(in.m_data))[inOffset+in.index(i1,i0)]) / 2.0;
887           }
888         }
889       }
890       if (in.getRank() == 4) {
891         int i0, i1, i2, i3;
892         int s0=in.getShape()[0];
893         int s1=in.getShape()[1];
894         int s2=in.getShape()[2];
895         int s3=in.getShape()[3];
896         for (i0=0; i0<s0; i0++) {
897           for (i1=0; i1<s1; i1++) {
898             for (i2=0; i2<s2; i2++) {
899               for (i3=0; i3<s3; i3++) {
900                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = ((*(in.m_data))[inOffset+in.index(i0,i1,i2,i3)] + (*(in.m_data))[inOffset+in.index(i2,i3,i0,i1)]) / 2.0;
901               }
902             }
903           }
904         }
905       }
906      }
907    
908      /**
909         \brief
910         computes a nonsymmetric matrix from your square matrix A: (A - transpose(A)) / 2
911    
912         \param in - Input - matrix
913         \param inOffset - Input - offset into in
914         \param ev - Output - The nonsymmetric matrix
915         \param inOffset - Input - offset into ev
916      */
917      static
918      inline
919      void
920      nonsymmetric(DataArrayView& in,
921                ValueType::size_type inOffset,
922                DataArrayView& ev,
923                ValueType::size_type evOffset)
924      {
925       if (in.getRank() == 2) {
926         int i0, i1;
927         int s0=in.getShape()[0];
928         int s1=in.getShape()[1];
929         for (i0=0; i0<s0; i0++) {
930           for (i1=0; i1<s1; i1++) {
931             (*(ev.m_data))[evOffset+ev.index(i0,i1)] = ((*(in.m_data))[inOffset+in.index(i0,i1)] - (*(in.m_data))[inOffset+in.index(i1,i0)]) / 2.0;
932           }
933         }
934       }
935       if (in.getRank() == 4) {
936         int i0, i1, i2, i3;
937         int s0=in.getShape()[0];
938         int s1=in.getShape()[1];
939         int s2=in.getShape()[2];
940         int s3=in.getShape()[3];
941         for (i0=0; i0<s0; i0++) {
942           for (i1=0; i1<s1; i1++) {
943             for (i2=0; i2<s2; i2++) {
944               for (i3=0; i3<s3; i3++) {
945                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = ((*(in.m_data))[inOffset+in.index(i0,i1,i2,i3)] - (*(in.m_data))[inOffset+in.index(i2,i3,i0,i1)]) / 2.0;
946               }
947             }
948           }
949         }
950       }
951      }
952    
953      /**
954         \brief
955         computes the trace of a matrix
956    
957         \param in - Input - matrix
958         \param inOffset - Input - offset into in
959         \param ev - Output - The nonsymmetric matrix
960         \param inOffset - Input - offset into ev
961      */
962      static
963      inline
964      void
965      matrixtrace(DataArrayView& in,
966                ValueType::size_type inOffset,
967                DataArrayView& ev,
968                ValueType::size_type evOffset,
969            int axis_offset)
970      {
971       if (in.getRank() == 2) {
972         int s0=in.getShape()[0]; // Python wrapper limits to square matrix
973         int i;
974         for (i=0; i<s0; i++) {
975           (*(ev.m_data))[evOffset+ev.index()] += (*(in.m_data))[inOffset+in.index(i,i)];
976         }
977       }
978       else if (in.getRank() == 3) {
979         if (axis_offset==0) {
980           int s0=in.getShape()[0];
981           int s2=in.getShape()[2];
982           int i0, i2;
983           for (i0=0; i0<s0; i0++) {
984             for (i2=0; i2<s2; i2++) {
985               (*(ev.m_data))[evOffset+ev.index(i2)] += (*(in.m_data))[inOffset+in.index(i0,i0,i2)];
986             }
987           }
988         }
989         else if (axis_offset==1) {
990           int s0=in.getShape()[0];
991           int s1=in.getShape()[1];
992           int i0, i1;
993           for (i0=0; i0<s0; i0++) {
994             for (i1=0; i1<s1; i1++) {
995               (*(ev.m_data))[evOffset+ev.index(i0)] += (*(in.m_data))[inOffset+in.index(i0,i1,i1)];
996             }
997           }
998         }
999       }
1000       else if (in.getRank() == 4) {
1001         if (axis_offset==0) {
1002           int s0=in.getShape()[0];
1003           int s2=in.getShape()[2];
1004           int s3=in.getShape()[3];
1005           int i0, i2, i3;
1006           for (i0=0; i0<s0; i0++) {
1007             for (i2=0; i2<s2; i2++) {
1008               for (i3=0; i3<s3; i3++) {
1009                 (*(ev.m_data))[evOffset+ev.index(i2,i3)] += (*(in.m_data))[inOffset+in.index(i0,i0,i2,i3)];
1010               }
1011             }
1012           }
1013         }
1014         else if (axis_offset==1) {
1015           int s0=in.getShape()[0];
1016           int s1=in.getShape()[1];
1017           int s3=in.getShape()[3];
1018           int i0, i1, i3;
1019           for (i0=0; i0<s0; i0++) {
1020             for (i1=0; i1<s1; i1++) {
1021               for (i3=0; i3<s3; i3++) {
1022                 (*(ev.m_data))[evOffset+ev.index(i0,i3)] += (*(in.m_data))[inOffset+in.index(i0,i1,i1,i3)];
1023               }
1024             }
1025           }
1026         }
1027         else if (axis_offset==2) {
1028           int s0=in.getShape()[0];
1029           int s1=in.getShape()[1];
1030           int s2=in.getShape()[2];
1031           int i0, i1, i2;
1032           for (i0=0; i0<s0; i0++) {
1033             for (i1=0; i1<s1; i1++) {
1034               for (i2=0; i2<s2; i2++) {
1035                 (*(ev.m_data))[evOffset+ev.index(i0,i1)] += (*(in.m_data))[inOffset+in.index(i0,i1,i2,i2)];
1036               }
1037             }
1038           }
1039         }
1040       }
1041      }
1042    
1043      /**
1044         \brief
1045         Transpose each data point of this Data object around the given axis.
1046    
1047         \param in - Input - matrix
1048         \param inOffset - Input - offset into in
1049         \param ev - Output - The nonsymmetric matrix
1050         \param inOffset - Input - offset into ev
1051      */
1052      static
1053      inline
1054      void
1055      transpose(DataArrayView& in,
1056                ValueType::size_type inOffset,
1057                DataArrayView& ev,
1058                ValueType::size_type evOffset,
1059            int axis_offset)
1060      {
1061       if (in.getRank() == 4) {
1062         int s0=ev.getShape()[0];
1063         int s1=ev.getShape()[1];
1064         int s2=ev.getShape()[2];
1065         int s3=ev.getShape()[3];
1066         int i0, i1, i2, i3;
1067         if (axis_offset==1) {
1068           for (i0=0; i0<s0; i0++) {
1069             for (i1=0; i1<s1; i1++) {
1070               for (i2=0; i2<s2; i2++) {
1071                 for (i3=0; i3<s3; i3++) {
1072                   (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i3,i0,i1,i2)];
1073                 }
1074               }
1075             }
1076           }
1077         }
1078         else if (axis_offset==2) {
1079           for (i0=0; i0<s0; i0++) {
1080             for (i1=0; i1<s1; i1++) {
1081               for (i2=0; i2<s2; i2++) {
1082                 for (i3=0; i3<s3; i3++) {
1083                   (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i2,i3,i0,i1)];
1084                 }
1085               }
1086             }
1087           }
1088         }
1089         else if (axis_offset==3) {
1090           for (i0=0; i0<s0; i0++) {
1091             for (i1=0; i1<s1; i1++) {
1092               for (i2=0; i2<s2; i2++) {
1093                 for (i3=0; i3<s3; i3++) {
1094                   (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i1,i2,i3,i0)];
1095                 }
1096               }
1097             }
1098           }
1099         }
1100         else {
1101           for (i0=0; i0<s0; i0++) {
1102             for (i1=0; i1<s1; i1++) {
1103               for (i2=0; i2<s2; i2++) {
1104                 for (i3=0; i3<s3; i3++) {
1105                   (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i0,i1,i2,i3)];
1106                 }
1107               }
1108             }
1109           }
1110         }
1111       }
1112       else if (in.getRank() == 3) {
1113         int s0=ev.getShape()[0];
1114         int s1=ev.getShape()[1];
1115         int s2=ev.getShape()[2];
1116         int i0, i1, i2;
1117         if (axis_offset==1) {
1118           for (i0=0; i0<s0; i0++) {
1119             for (i1=0; i1<s1; i1++) {
1120               for (i2=0; i2<s2; i2++) {
1121                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i2,i0,i1)];
1122               }
1123             }
1124           }
1125         }
1126         else if (axis_offset==2) {
1127           for (i0=0; i0<s0; i0++) {
1128             for (i1=0; i1<s1; i1++) {
1129               for (i2=0; i2<s2; i2++) {
1130                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i1,i2,i0)];
1131               }
1132             }
1133           }
1134         }
1135         else {
1136           // Copy the matrix unchanged
1137           for (i0=0; i0<s0; i0++) {
1138             for (i1=0; i1<s1; i1++) {
1139               for (i2=0; i2<s2; i2++) {
1140                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i0,i1,i2)];
1141               }
1142             }
1143           }
1144         }
1145       }
1146       else if (in.getRank() == 2) {
1147         int s0=ev.getShape()[0];
1148         int s1=ev.getShape()[1];
1149         int i0, i1;
1150         if (axis_offset==1) {
1151           for (i0=0; i0<s0; i0++) {
1152             for (i1=0; i1<s1; i1++) {
1153               (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i1,i0)];
1154             }
1155           }
1156         }
1157         else {
1158           for (i0=0; i0<s0; i0++) {
1159             for (i1=0; i1<s1; i1++) {
1160               (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i0,i1)];
1161             }
1162           }
1163         }
1164       }
1165       else if (in.getRank() == 1) {
1166         int s0=ev.getShape()[0];
1167         int i0;
1168         for (i0=0; i0<s0; i0++) {
1169           (*(ev.m_data))[evOffset+ev.index(i0)] = (*(in.m_data))[inOffset+in.index(i0)];
1170         }
1171       }
1172       else if (in.getRank() == 0) {
1173         (*(ev.m_data))[evOffset+ev.index()] = (*(in.m_data))[inOffset+in.index()];
1174       }
1175       else {
1176          throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
1177       }
1178      }
1179    
1180      /**
1181         \brief
1182       solves a local eigenvalue problem       solves a local eigenvalue problem
1183    
1184       \param in - Input - matrix       \param in - Input - matrix
Legend:
Removed from v.757  
changed lines
  Added in v.775
  ViewVC Help
Powered by ViewVC 1.1.26