Changeset 2575 in Sophya for trunk/SophyaLib/TArray/tarray.cc
- Timestamp:
- Jul 29, 2004, 2:31:16 PM (21 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/TArray/tarray.cc
r2569 r2575 482 482 //! Fill an array with a constant value \b x 483 483 template <class T> 484 TArray<T>& TArray<T>::Set T(T x)485 { 486 if (NbDimensions() < 1) 487 throw RangeCheckError("TArray<T>::Set T(T ) - Not Allocated Array ! ");484 TArray<T>& TArray<T>::SetCst(T x) 485 { 486 if (NbDimensions() < 1) 487 throw RangeCheckError("TArray<T>::SetCst(T ) - Not Allocated Array ! "); 488 488 T * pe; 489 489 sa_size_t j,k; … … 528 528 const T * pe; 529 529 T * per; 530 sa_size_t j,k ,kr;530 sa_size_t j,k; 531 531 if (smo && (IsPacked() > 0) && (res.IsPacked() > 0)) { // regularly spaced elements 532 532 sa_size_t maxx = totsize_; … … 538 538 int_4 ax,axr; 539 539 sa_size_t step, stepr; 540 sa_size_t gpas, nax r;541 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, nax r);542 for(j=0; j<nax r; j++) {540 sa_size_t gpas, naxa; 541 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa); 542 for(j=0; j<naxa; j++) { 543 543 pe = mNDBlock.Begin()+Offset(ax,j); 544 544 per = res.DataBlock().Begin()+res.Offset(axr,j); 545 for(k=0 , kr=0; k<gpas; k+=step, kr+=stepr) per[kr] = pe[k]+x;545 for(k=0; k<gpas; k+=step, pe+=step, per+=stepr) *per = *pe+x; 546 546 } 547 547 } … … 557 557 \param x : constant to subtract from the array elements 558 558 \param res : Output array containing the result (res=this+x or res=x-this). 559 \param fginv == true : Invert subtraction argument order ( *this = x-(*this))559 \param fginv == true : Invert subtraction argument order (res = x-(*this)) 560 560 */ 561 561 template <class T> … … 571 571 const T * pe; 572 572 T * per; 573 sa_size_t j,k ,kr;573 sa_size_t j,k; 574 574 if (smo && (IsPacked() > 0) && (res.IsPacked() > 0)) { // regularly spaced elements 575 575 sa_size_t maxx = totsize_; … … 584 584 int_4 ax,axr; 585 585 sa_size_t step, stepr; 586 sa_size_t gpas, nax r;587 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, nax r);588 for(j=0; j<nax r; j++) {586 sa_size_t gpas, naxa; 587 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa); 588 for(j=0; j<naxa; j++) { 589 589 pe = mNDBlock.Begin()+Offset(ax,j); 590 590 per = res.DataBlock().Begin()+res.Offset(axr,j); 591 591 if (!fginv) 592 for(k=0 , kr=0; k<gpas; k+=step, kr+=stepr) per[kr] = pe[k]-x;593 else594 for(k=0 , kr=0; k<gpas; k+=step, kr+=stepr) per[kr] = x-pe[k];592 for(k=0; k<gpas; k+=step, pe+=step, per+=stepr) *per = *pe-x; 593 else 594 for(k=0; k<gpas; k+=step, pe+=step, per+=stepr) *per = x-*pe; 595 595 } 596 596 } … … 619 619 const T * pe; 620 620 T * per; 621 sa_size_t j,k ,kr;621 sa_size_t j,k; 622 622 if (smo && (IsPacked() > 0) && (res.IsPacked() > 0)) { // regularly spaced elements 623 623 sa_size_t maxx = totsize_; … … 629 629 int_4 ax,axr; 630 630 sa_size_t step, stepr; 631 sa_size_t gpas, nax r;632 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, nax r);633 for(j=0; j<nax r; j++) {631 sa_size_t gpas, naxa; 632 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa); 633 for(j=0; j<naxa; j++) { 634 634 pe = mNDBlock.Begin()+Offset(ax,j); 635 635 per = res.DataBlock().Begin()+res.Offset(axr,j); 636 for(k=0 , kr=0; k<gpas; k+=step, kr+=stepr) per[kr] = pe[k]*x;636 for(k=0; k<gpas; k+=step, pe+=step, per+=stepr) *per = (*pe)*x; 637 637 } 638 638 } … … 648 648 \param x : Array elements are divied by x 649 649 \param res : Output array containing the result (res=(*this)/x or res=x/(*this)). 650 \param fginv == true : Invert the operationorder (res = x/(*this))650 \param fginv == true : Invert the division argument order (res = x/(*this)) 651 651 */ 652 652 template <class T> … … 677 677 int_4 ax,axr; 678 678 sa_size_t step, stepr; 679 sa_size_t gpas, nax r;680 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, nax r);681 for(j=0; j<nax r; j++) {679 sa_size_t gpas, naxa; 680 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa); 681 for(j=0; j<naxa; j++) { 682 682 pe = mNDBlock.Begin()+Offset(ax,j); 683 683 per = res.DataBlock().Begin()+res.Offset(axr,j); 684 684 if (!fginv) 685 for(k=0 , kr=0; k<gpas; k+=step, kr+=stepr) per[kr] = pe[k]/x;686 else687 for(k=0 , kr=0; k<gpas; k+=step, kr+=stepr) per[kr] = x/pe[k];685 for(k=0; k<gpas; k+=step, pe+=step, per+=stepr) *per = (*pe)/x; 686 else 687 for(k=0; k<gpas; k+=step, pe+=step, per+=stepr) *per = x/(*pe); 688 688 } 689 689 } … … 710 710 const T * pe; 711 711 T * per; 712 sa_size_t j,k ,kr;712 sa_size_t j,k; 713 713 if (smo && (IsPacked() > 0) && (res.IsPacked() > 0)) { // regularly spaced elements 714 714 sa_size_t maxx = totsize_; … … 720 720 int_4 ax,axr; 721 721 sa_size_t step, stepr; 722 sa_size_t gpas, nax r;723 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, nax r);724 for(j=0; j<nax r; j++) {722 sa_size_t gpas, naxa; 723 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa); 724 for(j=0; j<naxa; j++) { 725 725 pe = mNDBlock.Begin()+Offset(ax,j); 726 726 per = res.DataBlock().Begin()+res.Offset(axr,j); 727 for(k=0 , kr=0; k<gpas; k+=step, kr+=stepr) per[kr] = -pe[k];727 for(k=0; k<gpas; k+=step, pe+=step, per+=stepr) *per = -(*pe); 728 728 } 729 729 } … … 732 732 733 733 // >>>> Operations avec 2nd membre de type tableau 734 //! Add two TArrays 735 template <class T> 736 TArray<T>& TArray<T>::AddElt(const TArray<T>& a) 737 { 738 if (NbDimensions() < 1) 739 throw RangeCheckError("TArray<T>::AddElt(const TArray<T>& ) - Not Allocated Array ! "); 734 735 //! Two TArrays element by element addition 736 /*! 737 Perform element by element addition of the source array (this) and the \b a array 738 and store the result in \b res (res = *this+a). The source and argument arrays (this, a) 739 should have the same sizes. 740 If not initially allocated, the output array \b res is automatically 741 resized as a packed array with the same sizes as the source (this) array. 742 Returns a reference to the output array \b res. 743 \param a : Array to be added to the source array. 744 \param res : Output array containing the result (res=this+a). 745 */ 746 template <class T> 747 TArray<T>& TArray<T>::AddElt(const TArray<T>& a, TArray<T>& res) const 748 { 749 if (NbDimensions() < 1) 750 throw RangeCheckError("TArray<T>::AddElt(...) - Not allocated source array ! "); 751 bool smoa; 752 if (!CompareSizes(a, smoa)) 753 throw(SzMismatchError("TArray<T>::AddElt(...) SizeMismatch(this,a)")) ; 754 if (res.NbDimensions() < 1) res.SetSize(*this, true, false); 755 bool smor; 756 if (!CompareSizes(res, smor)) 757 throw(SzMismatchError("TArray<T>::AddElt(...) SizeMismatch(this,res) ")) ; 758 759 bool smora; 760 a.CompareSizes(res, smora); 761 762 bool smo = smoa && smor; // The three arrays have same memory organisation 763 764 const T * pe; 765 const T * pea; 766 T * per; 767 sa_size_t j,k; 768 if (smo && IsPacked() && a.IsPacked() && res.IsPacked() ) { // All packed arrays 769 sa_size_t maxx = totsize_; 770 pe = Data(); 771 pea = a.Data(); 772 per = res.Data(); 773 for(k=0; k<maxx; k++, pe++, pea++, per++) *per = *pe + *pea ; 774 } 775 else { // Non regular data spacing ... 776 int_4 ax,axa,axr; 777 sa_size_t step, stepa; 778 sa_size_t gpas, naxa; 779 sa_size_t stepr, stgpas; 780 if ( !smo && smora ) { // same mem-org for a,res , different from this 781 a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa); 782 a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa); 783 stgpas = stepa; 784 } 785 else { // same mem-org for all, or same (this,a) OR same(this,res) 786 GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa); 787 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa); 788 stgpas = step; 789 } 790 for(j=0; j<naxa; j++) { 791 pe = mNDBlock.Begin()+Offset(ax,j); 792 pea = a.DataBlock().Begin()+a.Offset(axa,j); 793 per = res.DataBlock().Begin()+res.Offset(axr,j); 794 for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = *pe + *pea ; 795 } 796 } 797 return(res); 798 } 799 800 //! Two TArrays element by element subtraction 801 /*! 802 Perform element by element subtraction of the source array (this) and the \b a array 803 and the store result in \b res (res = *this-a or res=a-(*this)). 804 The source and argument arrays (this, a) should have the same sizes. 805 If not initially allocated, the output array \b res is automatically 806 resized as a packed array with the same sizes as the source (this) array. 807 Returns a reference to the output array \b res. 808 \param a : Array to be added to the source array. 809 \param res : Output array containing the result (res=*this+x). 810 \param fginv == true : Invert subtraction argument order (res = a-(*this)) 811 */ 812 813 template <class T> 814 TArray<T>& TArray<T>::SubElt(const TArray<T>& a, TArray<T>& res, bool fginv) const 815 { 816 if (NbDimensions() < 1) 817 throw RangeCheckError("TArray<T>::SubElt(...) - Not allocated source array ! "); 818 bool smoa; 819 if (!CompareSizes(a, smoa)) 820 throw(SzMismatchError("TArray<T>::SubElt(...) SizeMismatch(this,a)")) ; 821 if (res.NbDimensions() < 1) res.SetSize(*this, true, false); 822 bool smor; 823 if (!CompareSizes(res, smor)) 824 throw(SzMismatchError("TArray<T>::SubElt(...) SizeMismatch(this,res) ")) ; 825 826 bool smora; 827 a.CompareSizes(res, smora); 828 829 bool smo = smoa && smor; // The three arrays have same memory organisation 830 831 const T * pe; 832 const T * pea; 833 T * per; 834 sa_size_t j,k; 835 if (smo && IsPacked() && a.IsPacked() && res.IsPacked() ) { // All packed arrays 836 sa_size_t maxx = totsize_; 837 pe = Data(); 838 pea = a.Data(); 839 per = res.Data(); 840 if (!fginv) 841 for(k=0; k<maxx; k++, pe++, pea++, per++) *per = *pe - *pea ; 842 else 843 for(k=0; k<maxx; k++, pe++, pea++, per++) *per = *pea - *pe ; 844 } 845 else { // Non regular data spacing ... 846 int_4 ax,axa,axr; 847 sa_size_t step, stepa; 848 sa_size_t gpas, naxa; 849 sa_size_t stepr, stgpas; 850 if ( !smo && smora ) { // same mem-org for a,res , different from this 851 a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa); 852 a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa); 853 stgpas = stepa; 854 } 855 else { // same mem-org for all, or same (this,a) OR same(this,res) 856 GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa); 857 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa); 858 stgpas = step; 859 } 860 for(j=0; j<naxa; j++) { 861 pe = mNDBlock.Begin()+Offset(ax,j); 862 pea = a.DataBlock().Begin()+a.Offset(axa,j); 863 per = res.DataBlock().Begin()+res.Offset(axr,j); 864 if (!fginv) 865 for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = *pe - *pea ; 866 else 867 for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = *pea - *pea ; 868 } 869 } 870 return(res); 871 } 872 873 874 //! Two TArrays element by element multiplication 875 /*! 876 Perform element by element multiplication of the source array (this) and the \b a array 877 and store the result in \b res (res = *this*a). The source and argument arrays (this, a) 878 should have the same sizes. 879 If not initially allocated, the output array \b res is automatically 880 resized as a packed array with the same sizes as the source (this) array. 881 Returns a reference to the output array \b res. 882 \param a : Array to be added to the source array. 883 \param res : Output array containing the result (res=(*this)*a). 884 */ 885 template <class T> 886 TArray<T>& TArray<T>::MulElt(const TArray<T>& a, TArray<T>& res) const 887 { 888 if (NbDimensions() < 1) 889 throw RangeCheckError("TArray<T>::MulElt(...) - Not allocated source array ! "); 890 bool smoa; 891 if (!CompareSizes(a, smoa)) 892 throw(SzMismatchError("TArray<T>::MulElt(...) SizeMismatch(this,a)")) ; 893 if (res.NbDimensions() < 1) res.SetSize(*this, true, false); 894 bool smor; 895 if (!CompareSizes(res, smor)) 896 throw(SzMismatchError("TArray<T>::MulElt(...) SizeMismatch(this,res) ")) ; 897 898 bool smora; 899 a.CompareSizes(res, smora); 900 901 bool smo = smoa && smor; // The three arrays have same memory organisation 902 903 const T * pe; 904 const T * pea; 905 T * per; 906 sa_size_t j,k; 907 if (smo && IsPacked() && a.IsPacked() && res.IsPacked() ) { // All packed arrays 908 sa_size_t maxx = totsize_; 909 pe = Data(); 910 pea = a.Data(); 911 per = res.Data(); 912 for(k=0; k<maxx; k++, pe++, pea++, per++ ) *per = (*pe) * (*pea) ; 913 } 914 else { // Non regular data spacing ... 915 int_4 ax,axa,axr; 916 sa_size_t step, stepa; 917 sa_size_t gpas, naxa; 918 sa_size_t stepr, stgpas; 919 if ( !smo && smora ) { // same mem-org for a,res , different from this 920 a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa); 921 a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa); 922 stgpas = stepa; 923 } 924 else { // same mem-org for all, or same (this,a) OR same(this,res) 925 GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa); 926 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa); 927 stgpas = step; 928 } 929 for(j=0; j<naxa; j++) { 930 pe = mNDBlock.Begin()+Offset(ax,j); 931 pea = a.DataBlock().Begin()+a.Offset(axa,j); 932 per = res.DataBlock().Begin()+res.Offset(axr,j); 933 for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) *per = (*pe) * (*pea); 934 } 935 } 936 return(res); 937 } 938 939 940 //! Two TArrays element by element division 941 /*! 942 Perform element by element division of the source array (this) and the \b a array 943 and store the result in \b res (res = *this/a). The source and argument arrays (this, a) 944 should have the same sizes. 945 If not initially allocated, the output array \b res is automatically 946 resized as a packed array with the same sizes as the source (this) array. 947 Returns a reference to the output array \b res. 948 \param a : Array to be added to the source array. 949 \param res : Output array containing the result (res=*this/a). 950 \param fginv == true : Inverts the division argument order (res = a/(*this)) 951 \param divzero == true : Result is set to zero (res(i)=0) if the operation's 952 second argument is equal to zero (a(i)/(*this)(i)==0) 953 */ 954 template <class T> 955 TArray<T>& TArray<T>::DivElt(const TArray<T>& a, TArray<T>& res, bool fginv, bool divzero) const 956 { 957 if (NbDimensions() < 1) 958 throw RangeCheckError("TArray<T>::DivElt(...) - Not allocated source array ! "); 959 bool smoa; 960 if (!CompareSizes(a, smoa)) 961 throw(SzMismatchError("TArray<T>::DivElt(...) SizeMismatch(this,a)")) ; 962 if (res.NbDimensions() < 1) res.SetSize(*this, true, false); 963 bool smor; 964 if (!CompareSizes(res, smor)) 965 throw(SzMismatchError("TArray<T>::DivElt(...) SizeMismatch(this,res) ")) ; 966 967 bool smora; 968 a.CompareSizes(res, smora); 969 970 bool smo = smoa && smor; // The three arrays have same memory organisation 971 972 const T * pe; 973 const T * pea; 974 T * per; 975 sa_size_t j,k; 976 if (smo && IsPacked() && a.IsPacked() && res.IsPacked() ) { // All packed arrays 977 sa_size_t maxx = totsize_; 978 pe = Data(); 979 pea = a.Data(); 980 per = res.Data(); 981 if(divzero) { 982 if (!fginv) 983 for(k=0; k<maxx; k++, pe++, pea++, per++ ) 984 if (*pea==(T)0) *per = (T)0; else *per = *pe / *pea ; 985 else 986 for(k=0; k<maxx; k++, pe++, pea++, per++ ) 987 if (*pe==(T)0) *per = (T)0; else *per = *pea / *pe ; 988 } 989 else { 990 if (!fginv) 991 for(k=0; k<maxx; k++, pe++, pea++, per++ ) *per = *pe / *pea ; 992 else 993 for(k=0; k<maxx; k++, pe++, pea++, per++ ) *per = *pea / *pe ; 994 } 995 } 996 else { // Non regular data spacing ... 997 int_4 ax,axa,axr; 998 sa_size_t step, stepa; 999 sa_size_t gpas, naxa; 1000 sa_size_t stepr, stgpas; 1001 if ( !smo && smora ) { // same mem-org for a,res , different from this 1002 a.GetOpeParams(*this, smo, axa, ax, stepa, step, gpas, naxa); 1003 a.GetOpeParams(res, smo, axa, axr, stepa, stepr, gpas, naxa); 1004 stgpas = stepa; 1005 } 1006 else { // same mem-org for all, or same (this,a) OR same(this,res) 1007 GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa); 1008 GetOpeParams(res, smo, ax, axr, step, stepr, gpas, naxa); 1009 stgpas = step; 1010 } 1011 // DBG cout << "DBG-A-DIVELT naxa=" << naxa << " gpas= " << gpas 1012 // << " step=" << step << " stepa=" << stepa << " stepr=" << stepr 1013 // << " ax= " << ax << " axa= " << axa << " axr= " << axr << endl; 1014 for(j=0; j<naxa; j++) { 1015 pe = mNDBlock.Begin()+Offset(ax,j); 1016 pea = a.DataBlock().Begin()+a.Offset(axa,j); 1017 per = res.DataBlock().Begin()+res.Offset(axr,j); 1018 if(divzero) { 1019 if (!fginv) 1020 for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) 1021 if (*pea==(T)0) *per = (T)0; else *per = *pe / *pea ; 1022 else 1023 for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) 1024 if (*pe==(T)0) *per = (T)0; else *per = *pea / *pe ; 1025 } 1026 else { 1027 if (!fginv) 1028 for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) 1029 *per = *pe / *pea ; 1030 else 1031 for(k=0; k<gpas; k+=stgpas, pe+=step, pea+=stepa, per+=stepr) 1032 *per = *pea / *pe ; 1033 } 1034 } 1035 } 1036 return(res); 1037 } 1038 1039 1040 //! Copy elements of \b a 1041 template <class T> 1042 TArray<T>& TArray<T>::CopyElt(const TArray<T>& a) 1043 { 1044 if (NbDimensions() < 1) 1045 throw RangeCheckError("TArray<T>::CopyElt(const TArray<T>& ) - Not Allocated Array ! "); 740 1046 bool smo; 741 1047 if (!CompareSizes(a, smo)) 742 throw(SzMismatchError("TArray<T>:: AddElt(const TArray<T>&) SizeMismatch")) ;1048 throw(SzMismatchError("TArray<T>::CopyElt(const TArray<T>&) SizeMismatch")) ; 743 1049 744 1050 T * pe; 745 1051 const T * pea; 746 sa_size_t j,k ,ka;747 if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) ) { // regularly spaced elements1052 sa_size_t j,k; 1053 if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) ) { // regularly spaced elements 748 1054 sa_size_t step = AvgStep(); 749 1055 sa_size_t stepa = a.AvgStep(); … … 751 1057 pe = Data(); 752 1058 pea = a.Data(); 753 for(k=0 , ka=0; k<maxx; k+=step, ka+=stepa ) pe[k] += pea[ka];1059 for(k=0; k<maxx; k+=step, pe+=step, pea+=stepa ) *pe = *pea ; 754 1060 } 755 1061 else { // Non regular data spacing ... … … 761 1067 pe = mNDBlock.Begin()+Offset(ax,j); 762 1068 pea = a.DataBlock().Begin()+a.Offset(axa,j); 763 for(k=0, ka=0; k<gpas; k+=step, ka+=stepa) pe[k] += pea[ka]; 764 } 765 } 766 return(*this); 767 } 768 769 //! Substract two TArrays 770 /*! 771 Substract two TArrays *this = *this-a 772 \param fginv == true : Perfoms the inverse subtraction (*this = a-(*this)) 773 */ 774 template <class T> 775 TArray<T>& TArray<T>::SubElt(const TArray<T>& a, bool fginv) 776 { 777 if (NbDimensions() < 1) 778 throw RangeCheckError("TArray<T>::SubElt(const TArray<T>& ) - Not Allocated Array ! "); 779 bool smo; 780 if (!CompareSizes(a, smo)) 781 throw(SzMismatchError("TArray<T>::SubElt(const TArray<T>&) SizeMismatch")) ; 782 783 T * pe; 784 const T * pea; 785 sa_size_t j,k,ka; 786 if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) ) { // regularly spaced elements 787 sa_size_t step = AvgStep(); 788 sa_size_t stepa = a.AvgStep(); 789 sa_size_t maxx = totsize_*step; 790 pe = Data(); 791 pea = a.Data(); 792 if (fginv) 793 for(k=0, ka=0; k<maxx; k+=step, ka+=stepa ) pe[k] = pea[ka]-pe[k] ; 794 else 795 for(k=0, ka=0; k<maxx; k+=step, ka+=stepa ) pe[k] -= pea[ka] ; 796 } 797 else { // Non regular data spacing ... 798 int_4 ax,axa; 799 sa_size_t step, stepa; 800 sa_size_t gpas, naxa; 801 GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa); 802 for(j=0; j<naxa; j++) { 803 pe = mNDBlock.Begin()+Offset(ax,j); 804 pea = a.DataBlock().Begin()+a.Offset(axa,j); 805 if (fginv) 806 for(k=0, ka=0; k<gpas; k+=step, ka+=stepa) pe[k] = pea[ka]-pe[k] ; 807 else 808 for(k=0, ka=0; k<gpas; k+=step, ka+=stepa) pe[k] -= pea[ka]; 809 } 810 } 811 return(*this); 812 } 813 814 815 //! Multiply two TArrays (elements by elements) 816 template <class T> 817 TArray<T>& TArray<T>::MulElt(const TArray<T>& a) 818 { 819 if (NbDimensions() < 1) 820 throw RangeCheckError("TArray<T>::MulElt(const TArray<T>& ) - Not Allocated Array ! "); 821 bool smo; 822 if (!CompareSizes(a, smo)) 823 throw(SzMismatchError("TArray<T>::MulElt(const TArray<T>&) SizeMismatch")) ; 824 825 T * pe; 826 const T * pea; 827 sa_size_t j,k,ka; 828 if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) ) { // regularly spaced elements 829 sa_size_t step = AvgStep(); 830 sa_size_t stepa = a.AvgStep(); 831 sa_size_t maxx = totsize_*step; 832 pe = Data(); 833 pea = a.Data(); 834 for(k=0, ka=0; k<maxx; k+=step, ka+=stepa ) pe[k] *= pea[ka] ; 835 } 836 else { // Non regular data spacing ... 837 int_4 ax,axa; 838 sa_size_t step, stepa; 839 sa_size_t gpas, naxa; 840 GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa); 841 for(j=0; j<naxa; j++) { 842 pe = mNDBlock.Begin()+Offset(axa,j); 843 pea = a.DataBlock().Begin()+a.Offset(axa,j); 844 for(k=0, ka=0; k<gpas; k+=step, ka+=stepa) pe[k] *= pea[ka]; 845 } 846 } 847 return(*this); 848 } 849 850 851 //! Divide two TArrays (elements by elements) 852 /*! 853 Divide two TArrays *this = (*this)/a 854 \param fginv == true : Perfoms the inverse division (*this = a/(*this)) 855 \param divzero == true : if a(i)==0. result is set to zero: (*this)(i)==0. 856 */ 857 template <class T> 858 TArray<T>& TArray<T>::DivElt(const TArray<T>& a, bool fginv, bool divzero) 859 { 860 if (NbDimensions() < 1) 861 throw RangeCheckError("TArray<T>::DivElt(const TArray<T>& ) - Not Allocated Array ! "); 862 bool smo; 863 if (!CompareSizes(a, smo)) 864 throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&) SizeMismatch")) ; 865 866 T * pe; 867 const T * pea; 868 sa_size_t j,k,ka; 869 if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) ) { // regularly spaced elements 870 sa_size_t step = AvgStep(); 871 sa_size_t stepa = a.AvgStep(); 872 sa_size_t maxx = totsize_*step; 873 pe = Data(); 874 pea = a.Data(); 875 if(divzero) { 876 if (fginv) 877 for(k=0, ka=0; k<maxx; k+=step, ka+=stepa ) 878 {if(pe[k]==(T)0) pe[k] = (T)0; else pe[k] = pea[ka]/pe[k];} 879 else 880 for(k=0, ka=0; k<maxx; k+=step, ka+=stepa ) 881 {if(pea[k]==(T)0) pe[k] = (T)0; else pe[k] /= pea[ka] ;} 882 } else { 883 if (fginv) 884 for(k=0, ka=0; k<maxx; k+=step, ka+=stepa ) pe[k] = pea[ka]/pe[k]; 885 else 886 for(k=0, ka=0; k<maxx; k+=step, ka+=stepa ) pe[k] /= pea[ka] ; 887 } 888 } 889 else { // Non regular data spacing ... 890 int_4 ax,axa; 891 sa_size_t step, stepa; 892 sa_size_t gpas, naxa; 893 GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa); 894 for(j=0; j<naxa; j++) { 895 pe = mNDBlock.Begin()+Offset(ax,j); 896 pea = a.DataBlock().Begin()+a.Offset(axa,j); 897 if(divzero) { 898 if (fginv) 899 for(k=0, ka=0; k<gpas; k+=step, ka+=stepa) 900 {if(pe[k]==(T)0) pe[k] = (T)0; else pe[k] = pea[ka]/pe[k];} 901 else 902 for(k=0, ka=0; k<gpas; k+=step, ka+=stepa) 903 {if(pea[k]==(T)0) pe[k] = (T)0; else pe[k] /= pea[ka];} 904 } else { 905 if (fginv) 906 for(k=0, ka=0; k<gpas; k+=step, ka+=stepa) pe[k] = pea[ka]/pe[k]; 907 else 908 for(k=0, ka=0; k<gpas; k+=step, ka+=stepa) pe[k] /= pea[ka]; 909 } 910 } 911 } 912 return(*this); 913 } 914 915 //! Copy elements of \b a 916 template <class T> 917 TArray<T>& TArray<T>::CopyElt(const TArray<T>& a) 918 { 919 if (NbDimensions() < 1) 920 throw RangeCheckError("TArray<T>::CopyElt(const TArray<T>& ) - Not Allocated Array ! "); 921 bool smo; 922 if (!CompareSizes(a, smo)) 923 throw(SzMismatchError("TArray<T>::CopyElt(const TArray<T>&) SizeMismatch")) ; 924 925 T * pe; 926 const T * pea; 927 sa_size_t j,k,ka; 928 if (smo && (AvgStep() > 0) && (a.AvgStep() > 0) ) { // regularly spaced elements 929 sa_size_t step = AvgStep(); 930 sa_size_t stepa = a.AvgStep(); 931 sa_size_t maxx = totsize_*step; 932 pe = Data(); 933 pea = a.Data(); 934 for(k=0, ka=0; k<maxx; k+=step, ka+=stepa ) pe[k] = pea[ka] ; 935 } 936 else { // Non regular data spacing ... 937 int_4 ax,axa; 938 sa_size_t step, stepa; 939 sa_size_t gpas, naxa; 940 GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa); 941 for(j=0; j<naxa; j++) { 942 pe = mNDBlock.Begin()+Offset(ax,j); 943 pea = a.DataBlock().Begin()+a.Offset(axa,j); 944 for(k=0, ka=0; k<gpas; k+=step, ka+=stepa) pe[k] = pea[ka]; 1069 for(k=0; k<gpas; k+=step, pe+=step, pea+=stepa) *pe = *pea; 945 1070 } 946 1071 } … … 986 1111 } 987 1112 1113 //! Return the the scalar product of the two arrays (Sum_k[(*this)(k)*a(k)]) 1114 template <class T> 1115 T TArray<T>::ScalarProduct(const TArray<T>& a) const 1116 { 1117 if (NbDimensions() < 1) 1118 throw RangeCheckError("TArray<T>::ScalarProduct(...) - Not allocated source array "); 1119 bool smo; 1120 if (!CompareSizes(a, smo)) 1121 throw(SzMismatchError("TArray<T>::ScalarProduct(...) SizeMismatch(this,a) ")) ; 1122 1123 T res = (T)(0); 1124 const T * pe; 1125 const T * pea; 1126 sa_size_t j,k; 1127 if (smo && (IsPacked() > 0) && (a.IsPacked() > 0)) { // regularly spaced elements 1128 sa_size_t maxx = totsize_; 1129 pe = Data(); 1130 pea = a.Data(); 1131 for(k=0; k<maxx; k++, pe++, pea++) res += (*pe)*(*pea); 1132 } 1133 else { // Non regular data spacing ... 1134 int_4 ax,axa; 1135 sa_size_t step, stepa; 1136 sa_size_t gpas, naxa; 1137 GetOpeParams(a, smo, ax, axa, step, stepa, gpas, naxa); 1138 for(j=0; j<naxa; j++) { 1139 pe = mNDBlock.Begin()+Offset(ax,j); 1140 pea = a.DataBlock().Begin()+a.Offset(axa,j); 1141 for(k=0; k<gpas; k+=step, pe+=step, pea+=stepa) res += (*pe)*(*pea); 1142 } 1143 } 1144 return(res); 1145 } 1146 988 1147 989 1148 // Somme et produit des elements 990 //! Sum all elements1149 //! Returns the sum of all array elements 991 1150 template <class T> 992 1151 T TArray<T>::Sum() const … … 1016 1175 } 1017 1176 1018 //! Multiplyall elements1177 //! Return the product of all elements 1019 1178 template <class T> 1020 1179 T TArray<T>::Product() const … … 1044 1203 } 1045 1204 1046 //! Returns the sum of all elements squared (Sum1205 //! Returns the sum of all array elements squared (Sum_k((*this)(k)*(*this)(k)). 1047 1206 template <class T> 1048 1207 T TArray<T>::SumX2() const
Note:
See TracChangeset
for help on using the changeset viewer.