001 // --- BEGIN LICENSE BLOCK ---
002 /*
003 * Copyright (c) 2009, Mikio L. Braun
004 * All rights reserved.
005 *
006 * Redistribution and use in source and binary forms, with or without
007 * modification, are permitted provided that the following conditions are
008 * met:
009 *
010 * * Redistributions of source code must retain the above copyright
011 * notice, this list of conditions and the following disclaimer.
012 *
013 * * Redistributions in binary form must reproduce the above
014 * copyright notice, this list of conditions and the following
015 * disclaimer in the documentation and/or other materials provided
016 * with the distribution.
017 *
018 * * Neither the name of the Technische Universit?t Berlin nor the
019 * names of its contributors may be used to endorse or promote
020 * products derived from this software without specific prior
021 * written permission.
022 *
023 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
024 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
025 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
026 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
027 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
028 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
029 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
030 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
031 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
032 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
033 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
034 */
035 // --- END LICENSE BLOCK ---
036
037 /*
038 * To change this template, choose Tools | Templates
039 * and open the template in the editor.
040 */
041 package org.jblas;
042
043 import org.jblas.exceptions.LapackException;
044
045 /**
046 * <p>Implementation of some Blas functions, mostly those which require linear runtime
047 * in the number of matrix elements. Because of the copying overhead when passing
048 * primitive arrays to native code, it doesn't make sense for these functions
049 * to be implemented in native code. The Java code is about as fast.</p>
050 *
051 * <p>The same conventions were used as in the native code, that is, for each array
052 * you also pass an index pointing to the starting index.</p>
053 *
054 * <p>These methods are mostly optimized for the case where the starting index is 0
055 * and the increment is 1.</p>
056 */
057 public class JavaBlas {
058
059 /** Exchange two vectors. */
060 public static void rswap(int n, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
061 if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
062 double z;
063 for (int i = 0; i < n; i++) {
064 z = dx[i];
065 dx[i] = dy[i];
066 dy[i] = z;
067 }
068 } else {
069 double z;
070 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; xi += incx, yi += incy, c++) {
071 z = dx[xi];
072 dx[xi] = dy[yi];
073 dy[yi] = z;
074 }
075 }
076 }
077
078 /** Copy dx to dy. */
079 public static void rcopy(int n, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
080 if (dxIdx < 0 || dxIdx + (n - 1) * incx >= dx.length) {
081 throw new LapackException("Java.raxpy", "Parameters for x aren't valid! (n = " + n + ", dx.length = " + dx.length + ", dxIdx = " + dxIdx + ", incx = " + incx + ")");
082 }
083 if (dyIdx < 0 || dyIdx + (n - 1) * incy >= dy.length) {
084 throw new LapackException("Java.raxpy", "Parameters for y aren't valid! (n = " + n + ", dy.length = " + dy.length + ", dyIdx = " + dyIdx + ", incy = " + incy + ")");
085 }
086 if (incx == 1 && incy == 1) {
087 System.arraycopy(dx, dxIdx, dy, dyIdx, n);
088 } else {
089 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; xi += incx, yi += incy, c++) {
090 dy[yi] = dx[xi];
091 }
092 }
093 }
094
095 /** Compute dy <- da * dx + dy. */
096 public static void raxpy(int n, double da, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
097 if (dxIdx < 0 || dxIdx + (n - 1) * incx >= dx.length) {
098 throw new LapackException("Java.raxpy", "Parameters for x aren't valid! (n = " + n + ", dx.length = " + dx.length + ", dxIdx = " + dxIdx + ", incx = " + incx + ")");
099 }
100 if (dyIdx < 0 || dyIdx + (n - 1) * incy >= dy.length) {
101 throw new LapackException("Java.raxpy", "Parameters for y aren't valid! (n = " + n + ", dy.length = " + dy.length + ", dyIdx = " + dyIdx + ", incy = " + incy + ")");
102 }
103
104 if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
105 if (da == 1.0) {
106 for (int i = 0; i < n; i++) {
107 dy[i] += dx[i];
108 }
109 } else {
110 for (int i = 0; i < n; i++) {
111 dy[i] += da * dx[i];
112 }
113 }
114 } else {
115 if (da == 1.0) {
116 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
117 dy[yi] += dx[xi];
118 }
119
120 } else {
121 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
122 dy[yi] += da * dx[xi];
123 }
124 }
125 }
126 }
127
128 /** Computes dz <- dx + dy */
129 public static void rzaxpy(int n, double[] dz, int dzIdx, int incz, double da, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
130 if (dxIdx == 0 && incx == 1 && dyIdx == 0 && incy == 1 && dzIdx == 0 && incz == 1) {
131 if (da == 1.0) {
132 for (int c = 0; c < n; c++)
133 dz[c] = dx[c] + dy[c];
134 } else {
135 for (int c = 0; c < n; c++)
136 dz[c] = da*dx[c] + dy[c];
137 }
138 } else {
139 if (da == 1.0) {
140 for (int c = 0, xi = dxIdx, yi = dyIdx, zi = dzIdx; c < n; c++, xi += incx, yi += incy, zi += incz) {
141 dz[zi] = dx[xi] + dy[yi];
142 }
143 } else {
144 for (int c = 0, xi = dxIdx, yi = dyIdx, zi = dzIdx; c < n; c++, xi += incx, yi += incy, zi += incz) {
145 dz[zi] = da*dx[xi] + dy[yi];
146 }
147 }
148 }
149 }
150
151 public static void rzgxpy(int n, double[] dz, double[] dx, double[] dy) {
152 for (int c = 0; c < n; c++)
153 dz[c] = dx[c] + dy[c];
154 }
155
156 /** Compute scalar product between dx and dy. */
157 public static double rdot(int n, double[] dx, int dxIdx, int incx, double[] dy, int dyIdx, int incy) {
158 double s = 0.0;
159 if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
160 for (int i = 0; i < n; i++)
161 s += dx[i] * dy[i];
162 }
163 else {
164 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
165 s += dx[xi] * dy[yi];
166 }
167 }
168 return s;
169 }
170 //BEGIN
171 // The code below has been automatically generated.
172 // DO NOT EDIT!
173
174 /** Exchange two vectors. */
175 public static void rswap(int n, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
176 if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
177 float z;
178 for (int i = 0; i < n; i++) {
179 z = dx[i];
180 dx[i] = dy[i];
181 dy[i] = z;
182 }
183 } else {
184 float z;
185 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; xi += incx, yi += incy, c++) {
186 z = dx[xi];
187 dx[xi] = dy[yi];
188 dy[yi] = z;
189 }
190 }
191 }
192
193 /** Copy dx to dy. */
194 public static void rcopy(int n, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
195 if (dxIdx < 0 || dxIdx + (n - 1) * incx >= dx.length) {
196 throw new LapackException("Java.raxpy", "Parameters for x aren't valid! (n = " + n + ", dx.length = " + dx.length + ", dxIdx = " + dxIdx + ", incx = " + incx + ")");
197 }
198 if (dyIdx < 0 || dyIdx + (n - 1) * incy >= dy.length) {
199 throw new LapackException("Java.raxpy", "Parameters for y aren't valid! (n = " + n + ", dy.length = " + dy.length + ", dyIdx = " + dyIdx + ", incy = " + incy + ")");
200 }
201 if (incx == 1 && incy == 1) {
202 System.arraycopy(dx, dxIdx, dy, dyIdx, n);
203 } else {
204 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; xi += incx, yi += incy, c++) {
205 dy[yi] = dx[xi];
206 }
207 }
208 }
209
210 /** Compute dy <- da * dx + dy. */
211 public static void raxpy(int n, float da, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
212 if (dxIdx < 0 || dxIdx + (n - 1) * incx >= dx.length) {
213 throw new LapackException("Java.raxpy", "Parameters for x aren't valid! (n = " + n + ", dx.length = " + dx.length + ", dxIdx = " + dxIdx + ", incx = " + incx + ")");
214 }
215 if (dyIdx < 0 || dyIdx + (n - 1) * incy >= dy.length) {
216 throw new LapackException("Java.raxpy", "Parameters for y aren't valid! (n = " + n + ", dy.length = " + dy.length + ", dyIdx = " + dyIdx + ", incy = " + incy + ")");
217 }
218
219 if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
220 if (da == 1.0f) {
221 for (int i = 0; i < n; i++) {
222 dy[i] += dx[i];
223 }
224 } else {
225 for (int i = 0; i < n; i++) {
226 dy[i] += da * dx[i];
227 }
228 }
229 } else {
230 if (da == 1.0f) {
231 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
232 dy[yi] += dx[xi];
233 }
234
235 } else {
236 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
237 dy[yi] += da * dx[xi];
238 }
239 }
240 }
241 }
242
243 /** Computes dz <- dx + dy */
244 public static void rzaxpy(int n, float[] dz, int dzIdx, int incz, float da, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
245 if (dxIdx == 0 && incx == 1 && dyIdx == 0 && incy == 1 && dzIdx == 0 && incz == 1) {
246 if (da == 1.0f) {
247 for (int c = 0; c < n; c++)
248 dz[c] = dx[c] + dy[c];
249 } else {
250 for (int c = 0; c < n; c++)
251 dz[c] = da*dx[c] + dy[c];
252 }
253 } else {
254 if (da == 1.0f) {
255 for (int c = 0, xi = dxIdx, yi = dyIdx, zi = dzIdx; c < n; c++, xi += incx, yi += incy, zi += incz) {
256 dz[zi] = dx[xi] + dy[yi];
257 }
258 } else {
259 for (int c = 0, xi = dxIdx, yi = dyIdx, zi = dzIdx; c < n; c++, xi += incx, yi += incy, zi += incz) {
260 dz[zi] = da*dx[xi] + dy[yi];
261 }
262 }
263 }
264 }
265
266 public static void rzgxpy(int n, float[] dz, float[] dx, float[] dy) {
267 for (int c = 0; c < n; c++)
268 dz[c] = dx[c] + dy[c];
269 }
270
271 /** Compute scalar product between dx and dy. */
272 public static float rdot(int n, float[] dx, int dxIdx, int incx, float[] dy, int dyIdx, int incy) {
273 float s = 0.0f;
274 if (incx == 1 && incy == 1 && dxIdx == 0 && dyIdx == 0) {
275 for (int i = 0; i < n; i++)
276 s += dx[i] * dy[i];
277 }
278 else {
279 for (int c = 0, xi = dxIdx, yi = dyIdx; c < n; c++, xi += incx, yi += incy) {
280 s += dx[xi] * dy[yi];
281 }
282 }
283 return s;
284 }
285 //END
286 }