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 }