001    package numal.highprecision;
002    
003    import numal.highprecision.*;
004    
005    /** A library of basic operations between matrices and vectrors.
006     * These methods make use of lower and upper indices; which allows
007     * vectors to be 0-index or 1-index (i.e. a[0...n-1] or a[1...n])
008     */
009    public class Basic extends Object {
010    
011        static final int BASE = 100;
012        
013        /** Compute a constant multiple (<i>x</i) of part of a vector <i>b</i>.
014         * Method originally from class numal.LinearAlgebra - Lau (2004) 
015         * section 1.2.A Real vector and matrix - Duplication p.6.
016         *
017         * @param l lower row bound index
018         * @param u upper row bound index
019         * @param shift index shifting parameter 
020         * @param a rectangular output matrix
021         * @param b rectangular input matrix
022         */
023        public static void dupvec(int l, int u, int shift, double a[], double b[]) {
024          for (; l<=u; l++) a[l]=b[l+shift];
025        }
026        
027        /** Replace a column sequence of elements of a rectangular matrix <i>a</i>
028         *  by a constant multiple (<i>x</i) of a column sequence of elements of a rectangular matrix <i>b</i>.
029         * Method originally from class numal.LinearAlgebra - Lau (2004) 
030         * section 1.3.C Real vector and matrix - Multiplication p.9.
031         *
032         * @param l lower row bound
033         * @param u upper row bound
034         * @param i column index of <i>a</i>
035         * @param j column index of <i>b</i>
036         * @param a rectangular matrix
037         * @param b rectangular matrix
038         * @param x multiplication factor
039         */
040        public static void mulcol(int l, int u, int i, int j, double a[][], double b[][], double x) {
041          for (; l<=u; l++) a[l][i]=b[l][j]*x;
042        }
043        
044        /** Compute the inner product of part of a vector <i>a</i> and part of a vector <i>b</i>
045         * A subset of the inner product
046         * can be computed by choosing suitable lower and upper bounds.
047         * Method originally from class numal.Basic - Lau (2004) 
048         * section 1.4.A Real vector vector products p.10.
049         *
050         * @param l lower bound
051         * @param u upper bound 
052         * @param shift index-shifting parameter of vector <i>b</i>
053         * @param a vector
054         * @param b vector
055         * @return the inner product of vector <i>a</i> and vector <i>b</i>
056         */
057        public static double vecvec(int l, int u, int shift, double a[], double b[]) {
058          int k;
059          double s;
060        
061          s=0.0;
062          for (k=l; k<=u; k++) s += a[k]*b[k+shift];
063          return (s);
064        }
065        
066        /** Compute the inner product of part of a row of a rectangular matrix <i>a</i>
067         * and the corresponding part of a vector <i>b</i>. A subset of the inner product
068         * can be computed by choosing suitable lower and upper column bounds.
069         * Method originally from class numal.Basic - Lau (2004) 
070         * section 1.4.B Real Vector Vector products p.11.
071         *
072         * @param l lower bound
073         * @param u upper bound 
074         * @param i row index of <i>a</i>
075         * @param a rectangular matrix
076         * @param b vector
077         * @return the inner product of row <i>i</i> of matrix <i>a</i> with vector <i>b</i>
078         */
079        public static double matvec(int l, int u, int i, double a[][], double b[]) {
080          int k;
081          double s;
082        
083          s=0.0;
084          for (k=l; k<=u; k++) s += a[i][k]*b[k];
085          return (s);
086        }
087        
088        /**
089         * Method originally from class numal.Basic - Lau (2004) 
090         * section 1.4.B Real Vector Vector products p.12.
091         * @param l lower bound
092         * @param u upper bound 
093         * @param i column index of <i>a</i>
094         * @param a rectangular matrix with column [l-u][j] defined
095         * @param b vector with entries [l-u] defined
096         */
097        public static double tamvec(int l, int u, int i, double a[][], double b[]) {
098          int k;
099          double s;
100        
101          s=0.0;
102          for (k=l; k<=u; k++) s += a[k][i]*b[k];
103          return (s);
104        }
105        
106        /** Get the inner product of part of a row of a rectangular matrix <i>a</i>
107         * and the corresponding part of a column of a rectangular matrix <i>b</i>.
108         *
109         * @param l lower row bound
110         * @param u upper row bound
111         * @param i row index of <i>a</i>
112         * @param j column index of <i>b</i>
113         * @param a rectangular matrix with row [i][l-u] defined
114         * @param b rectangular matrix with column [l-u][j] defined
115         * 
116         * @return inner product of <i>a</i> and <i>b</i>
117         */
118        public static double matmat(int l, int u, int i, int j, double a[][], double b[][]) {
119          int k;
120          double s;
121        
122          s=0.0;
123          for (k=l; k<=u; k++) s += a[i][k]*b[k][j];
124          return (s);
125        }
126        
127        /** Compute the inner product of part of a column of a rectangular matrix <i>a</i>
128         * and the corresponding part of a column of a rectangular matrix <i>b</i>.
129         * Method originally from class numal.Basic - Lau (2004) 
130         * section 1.4.E Real Vector Vector products p.13.
131         *
132         * @param l lower row bound
133         * @param u upper row bound
134         * @param i column index of <i>a</i>
135         * @param j column index of <i>b</i>
136         * @param a rectangular matrix with column [l-u][i] defined
137         * @param b rectangular matrix with column [l-u][j] defined
138         */
139        public static double tammat(int l, int u, int i, int j, double a[][], double b[][]) {
140          int k;
141          double s;
142        
143          s=0.0;
144          for (k=l; k<=u; k++) s += a[k][i]*b[k][j];
145          return (s);
146        }
147        
148        /** Compute the inner product of part of a row of a rectangular matrix <i>a</i> and the
149         * corresponding part of a row of a rectangular matrix <i>b</i>.
150         * Method originally from class numal.LinearAlgebra - Lau (2004) 
151         * section 1.4.F Real Vector Vector products p.13
152         *
153         * @param l lower column bound
154         * @param u upper column bound 
155         * @param i row index of <i>a</i>
156         * @param j row index of <i>b</i>
157         * @param a rectangular matrix
158         * @param b rectangular matrix
159         * @return inner product of row <i>i</i> of <i>a</i> and row <i>j</i> of <i>b</i>
160         */
161        public static double mattam(int l, int u, int i, int j, double a[][], double b[][]) {
162          int k;
163          double s;
164        
165          s=0.0;
166          for (k=l; k<=u; k++) s += a[i][k]*b[j][k];
167          return (s);
168        }
169        
170        /** Adds a constant multiple (<i>x</i>) of part of a column of a rectangular matrix
171         * <i>b</i> to part of a column of a rectangular matrix <i>a</i>
172         * Method originally from class numal.Basic - Lau (2004) 
173         * section 1.7.B Real vector and matrix - Elimination p.22.
174         *
175         * @param l lower row bound
176         * @param u upper row bound 
177         * @param i column index of <i>a</i>
178         * @param j column index of <i>b</i>
179         * @param a rectangular matrix
180         * @param b rectangular matrix
181         * @param x multiplication factor
182         */
183        public static void elmcol(int l, int u, int i, int j, double a[][], double b[][], double x) {
184            for (; l<=u; l++) a[l][i] += b[l][j]*x;
185        }
186        
187        /** Adds a constant multiple (<i>x</i>) of part of a row of a rectangular matrix
188         * <i>b</i> to part of a row of a rectangular matrix <i>a</i>
189         * Method originally from class numal.Basic - Lau (2004) 
190         * section 1.7.B Real vector and matrix - Elimination p.22.
191         *
192         * @param l lower row bound
193         * @param u upper row bound 
194         * @param i column index of <i>a</i>
195         * @param j column index of <i>b</i>
196         * @param a rectangular matrix
197         * @param b rectangular matrix
198         * @param x multiplication factor
199         */
200        public static void elmrow(int l, int u, int i, int j, double a[][], double b[][], double x) {
201          for (; l<=u; l++) a[i][l] += b[j][l]*x;
202        }
203        
204        /** Rotates two columns of a rectangular matrix using two parameters
205         *  <ul>
206         *  <li>a[k][i] = c * a[k][i] + s * a[k][j] (where k = l to u)</li>
207         *  <li>a[k][j] = c * a[k][j] - s * a[k][i]</li>
208         *  </ul>
209         * Method originally from class numal.Basic - Lau (2004) 
210         * section 1.9.A Real vector and matrix - Rotation p.30.
211         *
212         * @param l lower row bound
213         * @param u upper row bound 
214         * @param i column index of <i>a</i>
215         * @param j column index of <i>a</i>
216         * @param a rectangular matrix
217         * @param c cosine multiplication parameter
218         * @param s sine multiplication parameter
219         */
220        public static void rotcol(int l, int u, int i, int j, double a[][], double c, double s) {
221          double x, y;
222        
223          for (; l<=u; l++) {
224            x=a[l][i];
225            y=a[l][j];
226            a[l][i]=x*c+y*s;
227            a[l][j]=y*c-x*s;
228          }
229        }
230    
231    }