1 /*************************************************************************
2 Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
3
4 Redistribution and use in source and binary forms, with or without
5 modification, are permitted provided that the following conditions are
6 met:
7
8 - Redistributions of source code must retain the above copyright
9 notice, this list of conditions and the following disclaimer.
10
11 - Redistributions in binary form must reproduce the above copyright
12 notice, this list of conditions and the following disclaimer listed
13 in this license in the documentation and/or other materials
14 provided with the distribution.
15
16 - Neither the name of the copyright holders nor the names of its
17 contributors may be used to endorse or promote products derived from
18 this software without specific prior written permission.
19
20 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 *************************************************************************/
32
33 #include "alglib/blas.h"
34
vectornorm2(const ap::real_1d_array & x,int i1,int i2)35 double vectornorm2(const ap::real_1d_array& x, int i1, int i2)
36 {
37 double result;
38 int n;
39 int ix;
40 double absxi;
41 double scl;
42 double ssq;
43
44 n = i2-i1+1;
45 if( n<1 )
46 {
47 result = 0;
48 return result;
49 }
50 if( n==1 )
51 {
52 result = fabs(x(i1));
53 return result;
54 }
55 scl = 0;
56 ssq = 1;
57 for(ix = i1; ix <= i2; ix++)
58 {
59 if( x(ix)!=0 )
60 {
61 absxi = fabs(x(ix));
62 if( scl<absxi )
63 {
64 ssq = 1+ssq*ap::sqr(scl/absxi);
65 scl = absxi;
66 }
67 else
68 {
69 ssq = ssq+ap::sqr(absxi/scl);
70 }
71 }
72 }
73 result = scl*sqrt(ssq);
74 return result;
75 }
76
77
vectoridxabsmax(const ap::real_1d_array & x,int i1,int i2)78 int vectoridxabsmax(const ap::real_1d_array& x, int i1, int i2)
79 {
80 int result;
81 int i;
82
83 result = i1;
84 for(i = i1+1; i <= i2; i++)
85 {
86 if( fabs(x(i))>fabs(x(result)) )
87 {
88 result = i;
89 }
90 }
91 return result;
92 }
93
94
columnidxabsmax(const ap::real_2d_array & x,int i1,int i2,int j)95 int columnidxabsmax(const ap::real_2d_array& x, int i1, int i2, int j)
96 {
97 int result;
98 int i;
99
100 result = i1;
101 for(i = i1+1; i <= i2; i++)
102 {
103 if( fabs(x(i,j))>fabs(x(result,j)) )
104 {
105 result = i;
106 }
107 }
108 return result;
109 }
110
111
rowidxabsmax(const ap::real_2d_array & x,int j1,int j2,int i)112 int rowidxabsmax(const ap::real_2d_array& x, int j1, int j2, int i)
113 {
114 int result;
115 int j;
116
117 result = j1;
118 for(j = j1+1; j <= j2; j++)
119 {
120 if( fabs(x(i,j))>fabs(x(i,result)) )
121 {
122 result = j;
123 }
124 }
125 return result;
126 }
127
128
upperhessenberg1norm(const ap::real_2d_array & a,int i1,int i2,int j1,int j2,ap::real_1d_array & work)129 double upperhessenberg1norm(const ap::real_2d_array& a,
130 int i1,
131 int i2,
132 int j1,
133 int j2,
134 ap::real_1d_array& work)
135 {
136 double result;
137 int i;
138 int j;
139
140 ap::ap_error::make_assertion(i2-i1==j2-j1, "UpperHessenberg1Norm: I2-I1<>J2-J1!");
141 for(j = j1; j <= j2; j++)
142 {
143 work(j) = 0;
144 }
145 for(i = i1; i <= i2; i++)
146 {
147 for(j = ap::maxint(j1, j1+i-i1-1); j <= j2; j++)
148 {
149 work(j) = work(j)+fabs(a(i,j));
150 }
151 }
152 result = 0;
153 for(j = j1; j <= j2; j++)
154 {
155 result = ap::maxreal(result, work(j));
156 }
157 return result;
158 }
159
160
copymatrix(const ap::real_2d_array & a,int is1,int is2,int js1,int js2,ap::real_2d_array & b,int id1,int id2,int jd1,int jd2)161 void copymatrix(const ap::real_2d_array& a,
162 int is1,
163 int is2,
164 int js1,
165 int js2,
166 ap::real_2d_array& b,
167 int id1,
168 int id2,
169 int jd1,
170 int jd2)
171 {
172 int isrc;
173 int idst;
174
175 if( is1>is2||js1>js2 )
176 {
177 return;
178 }
179 ap::ap_error::make_assertion(is2-is1==id2-id1, "CopyMatrix: different sizes!");
180 ap::ap_error::make_assertion(js2-js1==jd2-jd1, "CopyMatrix: different sizes!");
181 for(isrc = is1; isrc <= is2; isrc++)
182 {
183 idst = isrc-is1+id1;
184 ap::vmove(&b(idst, jd1), &a(isrc, js1), ap::vlen(jd1,jd2));
185 }
186 }
187
188
inplacetranspose(ap::real_2d_array & a,int i1,int i2,int j1,int j2,ap::real_1d_array & work)189 void inplacetranspose(ap::real_2d_array& a,
190 int i1,
191 int i2,
192 int j1,
193 int j2,
194 ap::real_1d_array& work)
195 {
196 int i;
197 int j;
198 int ips;
199 int jps;
200 int l;
201
202 if( i1>i2||j1>j2 )
203 {
204 return;
205 }
206 ap::ap_error::make_assertion(i1-i2==j1-j2, "InplaceTranspose error: incorrect array size!");
207 for(i = i1; i <= i2-1; i++)
208 {
209 j = j1+i-i1;
210 ips = i+1;
211 jps = j1+ips-i1;
212 l = i2-i;
213 ap::vmove(work.getvector(1, l), a.getcolumn(j, ips, i2));
214 ap::vmove(a.getcolumn(j, ips, i2), a.getrow(i, jps, j2));
215 ap::vmove(&a(i, jps), &work(1), ap::vlen(jps,j2));
216 }
217 }
218
219
copyandtranspose(const ap::real_2d_array & a,int is1,int is2,int js1,int js2,ap::real_2d_array & b,int id1,int id2,int jd1,int jd2)220 void copyandtranspose(const ap::real_2d_array& a,
221 int is1,
222 int is2,
223 int js1,
224 int js2,
225 ap::real_2d_array& b,
226 int id1,
227 int id2,
228 int jd1,
229 int jd2)
230 {
231 int isrc;
232 int jdst;
233
234 if( is1>is2||js1>js2 )
235 {
236 return;
237 }
238 ap::ap_error::make_assertion(is2-is1==jd2-jd1, "CopyAndTranspose: different sizes!");
239 ap::ap_error::make_assertion(js2-js1==id2-id1, "CopyAndTranspose: different sizes!");
240 for(isrc = is1; isrc <= is2; isrc++)
241 {
242 jdst = isrc-is1+jd1;
243 ap::vmove(b.getcolumn(jdst, id1, id2), a.getrow(isrc, js1, js2));
244 }
245 }
246
247
matrixvectormultiply(const ap::real_2d_array & a,int i1,int i2,int j1,int j2,bool trans,const ap::real_1d_array & x,int ix1,int ix2,double alpha,ap::real_1d_array & y,int iy1,int iy2,double beta)248 void matrixvectormultiply(const ap::real_2d_array& a,
249 int i1,
250 int i2,
251 int j1,
252 int j2,
253 bool trans,
254 const ap::real_1d_array& x,
255 int ix1,
256 int ix2,
257 double alpha,
258 ap::real_1d_array& y,
259 int iy1,
260 int iy2,
261 double beta)
262 {
263 int i;
264 double v;
265
266 if( !trans )
267 {
268
269 //
270 // y := alpha*A*x + beta*y;
271 //
272 if( i1>i2||j1>j2 )
273 {
274 return;
275 }
276 ap::ap_error::make_assertion(j2-j1==ix2-ix1, "MatrixVectorMultiply: A and X dont match!");
277 ap::ap_error::make_assertion(i2-i1==iy2-iy1, "MatrixVectorMultiply: A and Y dont match!");
278
279 //
280 // beta*y
281 //
282 if( beta==0 )
283 {
284 for(i = iy1; i <= iy2; i++)
285 {
286 y(i) = 0;
287 }
288 }
289 else
290 {
291 ap::vmul(&y(iy1), ap::vlen(iy1,iy2), beta);
292 }
293
294 //
295 // alpha*A*x
296 //
297 for(i = i1; i <= i2; i++)
298 {
299 v = ap::vdotproduct(&a(i, j1), &x(ix1), ap::vlen(j1,j2));
300 y(iy1+i-i1) = y(iy1+i-i1)+alpha*v;
301 }
302 }
303 else
304 {
305
306 //
307 // y := alpha*A'*x + beta*y;
308 //
309 if( i1>i2||j1>j2 )
310 {
311 return;
312 }
313 ap::ap_error::make_assertion(i2-i1==ix2-ix1, "MatrixVectorMultiply: A and X dont match!");
314 ap::ap_error::make_assertion(j2-j1==iy2-iy1, "MatrixVectorMultiply: A and Y dont match!");
315
316 //
317 // beta*y
318 //
319 if( beta==0 )
320 {
321 for(i = iy1; i <= iy2; i++)
322 {
323 y(i) = 0;
324 }
325 }
326 else
327 {
328 ap::vmul(&y(iy1), ap::vlen(iy1,iy2), beta);
329 }
330
331 //
332 // alpha*A'*x
333 //
334 for(i = i1; i <= i2; i++)
335 {
336 v = alpha*x(ix1+i-i1);
337 ap::vadd(&y(iy1), &a(i, j1), ap::vlen(iy1,iy2), v);
338 }
339 }
340 }
341
342
pythag2(double x,double y)343 double pythag2(double x, double y)
344 {
345 double result;
346 double w;
347 double xabs;
348 double yabs;
349 double z;
350
351 xabs = fabs(x);
352 yabs = fabs(y);
353 w = ap::maxreal(xabs, yabs);
354 z = ap::minreal(xabs, yabs);
355 if( z==0 )
356 {
357 result = w;
358 }
359 else
360 {
361 result = w*sqrt(1+ap::sqr(z/w));
362 }
363 return result;
364 }
365
366
matrixmatrixmultiply(const ap::real_2d_array & a,int ai1,int ai2,int aj1,int aj2,bool transa,const ap::real_2d_array & b,int bi1,int bi2,int bj1,int bj2,bool transb,double alpha,ap::real_2d_array & c,int ci1,int ci2,int cj1,int cj2,double beta,ap::real_1d_array & work)367 void matrixmatrixmultiply(const ap::real_2d_array& a,
368 int ai1,
369 int ai2,
370 int aj1,
371 int aj2,
372 bool transa,
373 const ap::real_2d_array& b,
374 int bi1,
375 int bi2,
376 int bj1,
377 int bj2,
378 bool transb,
379 double alpha,
380 ap::real_2d_array& c,
381 int ci1,
382 int ci2,
383 int cj1,
384 int cj2,
385 double beta,
386 ap::real_1d_array& work)
387 {
388 int arows;
389 int acols;
390 int brows;
391 int bcols;
392 int crows;
393 int i;
394 int j;
395 int k = 0; // Eliminate compiler warning.
396 int l;
397 int r;
398 double v;
399
400
401 //
402 // Setup
403 //
404 if( !transa )
405 {
406 arows = ai2-ai1+1;
407 acols = aj2-aj1+1;
408 }
409 else
410 {
411 arows = aj2-aj1+1;
412 acols = ai2-ai1+1;
413 }
414 if( !transb )
415 {
416 brows = bi2-bi1+1;
417 bcols = bj2-bj1+1;
418 }
419 else
420 {
421 brows = bj2-bj1+1;
422 bcols = bi2-bi1+1;
423 }
424 ap::ap_error::make_assertion(acols==brows, "MatrixMatrixMultiply: incorrect matrix sizes!");
425 if( arows<=0||acols<=0||brows<=0||bcols<=0 )
426 {
427 return;
428 }
429 crows = arows;
430
431 //
432 // Test WORK
433 //
434 i = ap::maxint(arows, acols);
435 i = ap::maxint(brows, i);
436 i = ap::maxint(i, bcols);
437 work(1) = 0;
438 work(i) = 0;
439
440 //
441 // Prepare C
442 //
443 if( beta==0 )
444 {
445 for(i = ci1; i <= ci2; i++)
446 {
447 for(j = cj1; j <= cj2; j++)
448 {
449 c(i,j) = 0;
450 }
451 }
452 }
453 else
454 {
455 for(i = ci1; i <= ci2; i++)
456 {
457 ap::vmul(&c(i, cj1), ap::vlen(cj1,cj2), beta);
458 }
459 }
460
461 //
462 // A*B
463 //
464 if( !transa&&!transb )
465 {
466 for(l = ai1; l <= ai2; l++)
467 {
468 for(r = bi1; r <= bi2; r++)
469 {
470 v = alpha*a(l,aj1+r-bi1);
471 k = ci1+l-ai1;
472 ap::vadd(&c(k, cj1), &b(r, bj1), ap::vlen(cj1,cj2), v);
473 }
474 }
475 return;
476 }
477
478 //
479 // A*B'
480 //
481 if( !transa&&transb )
482 {
483 if( arows*acols<brows*bcols )
484 {
485 for(r = bi1; r <= bi2; r++)
486 {
487 for(l = ai1; l <= ai2; l++)
488 {
489 v = ap::vdotproduct(&a(l, aj1), &b(r, bj1), ap::vlen(aj1,aj2));
490 c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*v;
491 }
492 }
493 return;
494 }
495 else
496 {
497 for(l = ai1; l <= ai2; l++)
498 {
499 for(r = bi1; r <= bi2; r++)
500 {
501 v = ap::vdotproduct(&a(l, aj1), &b(r, bj1), ap::vlen(aj1,aj2));
502 c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*v;
503 }
504 }
505 return;
506 }
507 }
508
509 //
510 // A'*B
511 //
512 if( transa&&!transb )
513 {
514 for(l = aj1; l <= aj2; l++)
515 {
516 for(r = bi1; r <= bi2; r++)
517 {
518 v = alpha*a(ai1+r-bi1,l);
519 k = ci1+l-aj1;
520 ap::vadd(&c(k, cj1), &b(r, bj1), ap::vlen(cj1,cj2), v);
521 }
522 }
523 return;
524 }
525
526 //
527 // A'*B'
528 //
529 if( transa&&transb )
530 {
531 if( arows*acols<brows*bcols )
532 {
533 for(r = bi1; r <= bi2; r++)
534 {
535 for(i = 1; i <= crows; i++)
536 {
537 work(i) = 0.0;
538 }
539 for(l = ai1; l <= ai2; l++)
540 {
541 v = alpha*b(r,bj1+l-ai1);
542 k = cj1+r-bi1;
543 ap::vadd(&work(1), &a(l, aj1), ap::vlen(1,crows), v);
544 }
545 ap::vadd(c.getcolumn(k, ci1, ci2), work.getvector(1, crows));
546 }
547 return;
548 }
549 else
550 {
551 for(l = aj1; l <= aj2; l++)
552 {
553 k = ai2-ai1+1;
554 ap::vmove(work.getvector(1, k), a.getcolumn(l, ai1, ai2));
555 for(r = bi1; r <= bi2; r++)
556 {
557 v = ap::vdotproduct(&work(1), &b(r, bj1), ap::vlen(1,k));
558 c(ci1+l-aj1,cj1+r-bi1) = c(ci1+l-aj1,cj1+r-bi1)+alpha*v;
559 }
560 }
561 return;
562 }
563 }
564 }
565
566
567
568