OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [ada/] [s-gearop.adb] - Blame information for rev 706

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 706 jeremybenn
------------------------------------------------------------------------------
2
--                                                                          --
3
--                         GNAT RUN-TIME COMPONENTS                         --
4
--                                                                          --
5
--       S Y S T E M . G E N E R I C _ A R R A Y _ O P E R A T I O N S      --
6
--                                                                          --
7
--                                 B o d y                                  --
8
--                                                                          --
9
--         Copyright (C) 2006-2012, Free Software Foundation, Inc.          --
10
--                                                                          --
11
-- GNAT is free software;  you can  redistribute it  and/or modify it under --
12
-- terms of the  GNU General Public License as published  by the Free Soft- --
13
-- ware  Foundation;  either version 3,  or (at your option) any later ver- --
14
-- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
15
-- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
16
-- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
17
--                                                                          --
18
-- As a special exception under Section 7 of GPL version 3, you are granted --
19
-- additional permissions described in the GCC Runtime Library Exception,   --
20
-- version 3.1, as published by the Free Software Foundation.               --
21
--                                                                          --
22
-- You should have received a copy of the GNU General Public License and    --
23
-- a copy of the GCC Runtime Library Exception along with this program;     --
24
-- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
25
-- <http://www.gnu.org/licenses/>.                                          --
26
--                                                                          --
27
-- GNAT was originally developed  by the GNAT team at  New York University. --
28
-- Extensive contributions were provided by Ada Core Technologies Inc.      --
29
--                                                                          --
30
------------------------------------------------------------------------------
31
 
32
with Ada.Numerics; use Ada.Numerics;
33
 
34
package body System.Generic_Array_Operations is
35
 
36
   --  The local function Check_Unit_Last computes the index of the last
37
   --  element returned by Unit_Vector or Unit_Matrix. A separate function is
38
   --  needed to allow raising Constraint_Error before declaring the function
39
   --  result variable. The result variable needs to be declared first, to
40
   --  allow front-end inlining.
41
 
42
   function Check_Unit_Last
43
     (Index : Integer;
44
      Order : Positive;
45
      First : Integer) return Integer;
46
   pragma Inline_Always (Check_Unit_Last);
47
 
48
   --------------
49
   -- Diagonal --
50
   --------------
51
 
52
   function Diagonal (A : Matrix) return Vector is
53
      N : constant Natural := Natural'Min (A'Length (1), A'Length (2));
54
      R : Vector (A'First (1) .. A'First (1) + N - 1);
55
 
56
   begin
57
      for J in 0 .. N - 1 loop
58
         R (R'First + J) := A (A'First (1) + J, A'First (2) + J);
59
      end loop;
60
 
61
      return R;
62
   end Diagonal;
63
 
64
   --------------------------
65
   -- Square_Matrix_Length --
66
   --------------------------
67
 
68
   function Square_Matrix_Length (A : Matrix) return Natural is
69
   begin
70
      if A'Length (1) /= A'Length (2) then
71
         raise Constraint_Error with "matrix is not square";
72
      end if;
73
 
74
      return A'Length (1);
75
   end Square_Matrix_Length;
76
 
77
   ---------------------
78
   -- Check_Unit_Last --
79
   ---------------------
80
 
81
   function Check_Unit_Last
82
      (Index : Integer;
83
       Order : Positive;
84
       First : Integer) return Integer
85
   is
86
   begin
87
      --  Order the tests carefully to avoid overflow
88
 
89
      if Index < First
90
        or else First > Integer'Last - Order + 1
91
        or else Index > First + (Order - 1)
92
      then
93
         raise Constraint_Error;
94
      end if;
95
 
96
      return First + (Order - 1);
97
   end Check_Unit_Last;
98
 
99
   ---------------------
100
   -- Back_Substitute --
101
   ---------------------
102
 
103
   procedure Back_Substitute (M, N : in out Matrix) is
104
      pragma Assert (M'First (1) = N'First (1)
105
                       and then
106
                     M'Last  (1) = N'Last (1));
107
 
108
      procedure Sub_Row
109
        (M      : in out Matrix;
110
         Target : Integer;
111
         Source : Integer;
112
         Factor : Scalar);
113
      --  Elementary row operation that subtracts Factor * M (Source, <>) from
114
      --  M (Target, <>)
115
 
116
      procedure Sub_Row
117
        (M      : in out Matrix;
118
         Target : Integer;
119
         Source : Integer;
120
         Factor : Scalar)
121
      is
122
      begin
123
         for J in M'Range (2) loop
124
            M (Target, J) := M (Target, J) - Factor * M (Source, J);
125
         end loop;
126
      end Sub_Row;
127
 
128
      --  Local declarations
129
 
130
      Max_Col : Integer := M'Last (2);
131
 
132
   --  Start of processing for Back_Substitute
133
 
134
   begin
135
      Do_Rows : for Row in reverse M'Range (1) loop
136
         Find_Non_Zero : for Col in reverse M'First (2) .. Max_Col loop
137
            if Is_Non_Zero (M (Row, Col)) then
138
 
139
               --  Found first non-zero element, so subtract a multiple of this
140
               --  element  from all higher rows, to reduce all other elements
141
               --  in this column to zero.
142
 
143
               declare
144
                  --  We can't use a for loop, as we'd need to iterate to
145
                  --  Row - 1, but that expression will overflow if M'First
146
                  --  equals Integer'First, which is true for aggregates
147
                  --  without explicit bounds..
148
 
149
                  J : Integer := M'First (1);
150
 
151
               begin
152
                  while J < Row loop
153
                     Sub_Row (N, J, Row, (M (J, Col) / M (Row, Col)));
154
                     Sub_Row (M, J, Row, (M (J, Col) / M (Row, Col)));
155
                     J := J + 1;
156
                  end loop;
157
               end;
158
 
159
               --  Avoid potential overflow in the subtraction below
160
 
161
               exit Do_Rows when Col = M'First (2);
162
 
163
               Max_Col := Col - 1;
164
 
165
               exit Find_Non_Zero;
166
            end if;
167
         end loop Find_Non_Zero;
168
      end loop Do_Rows;
169
   end Back_Substitute;
170
 
171
   -----------------------
172
   -- Forward_Eliminate --
173
   -----------------------
174
 
175
   procedure Forward_Eliminate
176
     (M   : in out Matrix;
177
      N   : in out Matrix;
178
      Det : out Scalar)
179
   is
180
      pragma Assert (M'First (1) = N'First (1)
181
                       and then
182
                     M'Last  (1) = N'Last (1));
183
 
184
      --  The following are variations of the elementary matrix row operations:
185
      --  row switching, row multiplication and row addition. Because in this
186
      --  algorithm the addition factor is always a negated value, we chose to
187
      --  use  row subtraction instead. Similarly, instead of multiplying by
188
      --  a reciprocal, we divide.
189
 
190
      procedure Sub_Row
191
        (M      : in out Matrix;
192
         Target : Integer;
193
         Source : Integer;
194
         Factor : Scalar);
195
      --  Subtrace Factor * M (Source, <>) from M (Target, <>)
196
 
197
      procedure Divide_Row
198
        (M, N  : in out Matrix;
199
         Row   : Integer;
200
         Scale : Scalar);
201
      --  Divide M (Row) and N (Row) by Scale, and update Det
202
 
203
      procedure Switch_Row
204
        (M, N  : in out Matrix;
205
         Row_1 : Integer;
206
         Row_2 : Integer);
207
      --  Exchange M (Row_1) and N (Row_1) with M (Row_2) and N (Row_2),
208
      --  negating Det in the process.
209
 
210
      -------------
211
      -- Sub_Row --
212
      -------------
213
 
214
      procedure Sub_Row
215
        (M      : in out Matrix;
216
         Target : Integer;
217
         Source : Integer;
218
         Factor : Scalar)
219
      is
220
      begin
221
         for J in M'Range (2) loop
222
            M (Target, J) := M (Target, J) - Factor * M (Source, J);
223
         end loop;
224
      end Sub_Row;
225
 
226
      ----------------
227
      -- Divide_Row --
228
      ----------------
229
 
230
      procedure Divide_Row
231
        (M, N  : in out Matrix;
232
         Row   : Integer;
233
         Scale : Scalar)
234
      is
235
      begin
236
         Det := Det * Scale;
237
 
238
         for J in M'Range (2) loop
239
            M (Row, J) := M (Row, J) / Scale;
240
         end loop;
241
 
242
         for J in N'Range (2) loop
243
            N (Row - M'First (1) + N'First (1), J) :=
244
              N (Row - M'First (1) + N'First (1), J) / Scale;
245
         end loop;
246
      end Divide_Row;
247
 
248
      ----------------
249
      -- Switch_Row --
250
      ----------------
251
 
252
      procedure Switch_Row
253
        (M, N  : in out Matrix;
254
         Row_1 : Integer;
255
         Row_2 : Integer)
256
      is
257
         procedure Swap (X, Y : in out Scalar);
258
         --  Exchange the values of X and Y
259
 
260
         procedure Swap (X, Y : in out Scalar) is
261
            T : constant Scalar := X;
262
         begin
263
            X := Y;
264
            Y := T;
265
         end Swap;
266
 
267
      --  Start of processing for Switch_Row
268
 
269
      begin
270
         if Row_1 /= Row_2 then
271
            Det := Zero - Det;
272
 
273
            for J in M'Range (2) loop
274
               Swap (M (Row_1, J), M (Row_2, J));
275
            end loop;
276
 
277
            for J in N'Range (2) loop
278
               Swap (N (Row_1 - M'First (1) + N'First (1), J),
279
                     N (Row_2 - M'First (1) + N'First (1), J));
280
            end loop;
281
         end if;
282
      end Switch_Row;
283
 
284
      --  Local declarations
285
 
286
      Row : Integer := M'First (1);
287
 
288
   --  Start of processing for Forward_Eliminate
289
 
290
   begin
291
      Det := One;
292
 
293
      for J in M'Range (2) loop
294
         declare
295
            Max_Row : Integer := Row;
296
            Max_Abs : Real'Base := 0.0;
297
 
298
         begin
299
            --  Find best pivot in column J, starting in row Row
300
 
301
            for K in Row .. M'Last (1) loop
302
               declare
303
                  New_Abs : constant Real'Base := abs M (K, J);
304
               begin
305
                  if Max_Abs < New_Abs then
306
                     Max_Abs := New_Abs;
307
                     Max_Row := K;
308
                  end if;
309
               end;
310
            end loop;
311
 
312
            if Max_Abs > 0.0 then
313
               Switch_Row (M, N, Row, Max_Row);
314
 
315
               --  The temporaries below are necessary to force a copy of the
316
               --  value and avoid improper aliasing.
317
 
318
               declare
319
                  Scale : constant Scalar := M (Row, J);
320
               begin
321
                  Divide_Row (M, N, Row, Scale);
322
               end;
323
 
324
               for U in Row + 1 .. M'Last (1) loop
325
                  declare
326
                     Factor : constant Scalar := M (U, J);
327
                  begin
328
                     Sub_Row (N, U, Row, Factor);
329
                     Sub_Row (M, U, Row, Factor);
330
                  end;
331
               end loop;
332
 
333
               exit when Row >= M'Last (1);
334
 
335
               Row := Row + 1;
336
 
337
            else
338
               --  Set zero (note that we do not have literals)
339
 
340
               Det := Zero;
341
            end if;
342
         end;
343
      end loop;
344
   end Forward_Eliminate;
345
 
346
   -------------------
347
   -- Inner_Product --
348
   -------------------
349
 
350
   function Inner_Product
351
     (Left  : Left_Vector;
352
      Right : Right_Vector) return  Result_Scalar
353
   is
354
      R : Result_Scalar := Zero;
355
 
356
   begin
357
      if Left'Length /= Right'Length then
358
         raise Constraint_Error with
359
            "vectors are of different length in inner product";
360
      end if;
361
 
362
      for J in Left'Range loop
363
         R := R + Left (J) * Right (J - Left'First + Right'First);
364
      end loop;
365
 
366
      return R;
367
   end Inner_Product;
368
 
369
   -------------
370
   -- L2_Norm --
371
   -------------
372
 
373
   function L2_Norm (X : X_Vector) return Result_Real'Base is
374
      Sum : Result_Real'Base := 0.0;
375
 
376
   begin
377
      for J in X'Range loop
378
         Sum := Sum + Result_Real'Base (abs X (J))**2;
379
      end loop;
380
 
381
      return Sqrt (Sum);
382
   end L2_Norm;
383
 
384
   ----------------------------------
385
   -- Matrix_Elementwise_Operation --
386
   ----------------------------------
387
 
388
   function Matrix_Elementwise_Operation (X : X_Matrix) return Result_Matrix is
389
      R : Result_Matrix (X'Range (1), X'Range (2));
390
 
391
   begin
392
      for J in R'Range (1) loop
393
         for K in R'Range (2) loop
394
            R (J, K) := Operation (X (J, K));
395
         end loop;
396
      end loop;
397
 
398
      return R;
399
   end Matrix_Elementwise_Operation;
400
 
401
   ----------------------------------
402
   -- Vector_Elementwise_Operation --
403
   ----------------------------------
404
 
405
   function Vector_Elementwise_Operation (X : X_Vector) return Result_Vector is
406
      R : Result_Vector (X'Range);
407
 
408
   begin
409
      for J in R'Range loop
410
         R (J) := Operation (X (J));
411
      end loop;
412
 
413
      return R;
414
   end Vector_Elementwise_Operation;
415
 
416
   -----------------------------------------
417
   -- Matrix_Matrix_Elementwise_Operation --
418
   -----------------------------------------
419
 
420
   function Matrix_Matrix_Elementwise_Operation
421
     (Left  : Left_Matrix;
422
      Right : Right_Matrix) return Result_Matrix
423
   is
424
      R : Result_Matrix (Left'Range (1), Left'Range (2));
425
 
426
   begin
427
      if Left'Length (1) /= Right'Length (1)
428
           or else
429
         Left'Length (2) /= Right'Length (2)
430
      then
431
         raise Constraint_Error with
432
           "matrices are of different dimension in elementwise operation";
433
      end if;
434
 
435
      for J in R'Range (1) loop
436
         for K in R'Range (2) loop
437
            R (J, K) :=
438
              Operation
439
                (Left (J, K),
440
                 Right
441
                   (J - R'First (1) + Right'First (1),
442
                    K - R'First (2) + Right'First (2)));
443
         end loop;
444
      end loop;
445
 
446
      return R;
447
   end Matrix_Matrix_Elementwise_Operation;
448
 
449
   ------------------------------------------------
450
   -- Matrix_Matrix_Scalar_Elementwise_Operation --
451
   ------------------------------------------------
452
 
453
   function Matrix_Matrix_Scalar_Elementwise_Operation
454
     (X    : X_Matrix;
455
      Y    : Y_Matrix;
456
      Z    : Z_Scalar) return Result_Matrix
457
   is
458
      R : Result_Matrix (X'Range (1), X'Range (2));
459
 
460
   begin
461
      if X'Length (1) /= Y'Length (1)
462
           or else
463
         X'Length (2) /= Y'Length (2)
464
      then
465
         raise Constraint_Error with
466
           "matrices are of different dimension in elementwise operation";
467
      end if;
468
 
469
      for J in R'Range (1) loop
470
         for K in R'Range (2) loop
471
            R (J, K) :=
472
              Operation
473
                (X (J, K),
474
                 Y (J - R'First (1) + Y'First (1),
475
                    K - R'First (2) + Y'First (2)),
476
                 Z);
477
         end loop;
478
      end loop;
479
 
480
      return R;
481
   end Matrix_Matrix_Scalar_Elementwise_Operation;
482
 
483
   -----------------------------------------
484
   -- Vector_Vector_Elementwise_Operation --
485
   -----------------------------------------
486
 
487
   function Vector_Vector_Elementwise_Operation
488
     (Left  : Left_Vector;
489
      Right : Right_Vector) return Result_Vector
490
   is
491
      R : Result_Vector (Left'Range);
492
 
493
   begin
494
      if Left'Length /= Right'Length then
495
         raise Constraint_Error with
496
           "vectors are of different length in elementwise operation";
497
      end if;
498
 
499
      for J in R'Range loop
500
         R (J) := Operation (Left (J), Right (J - R'First + Right'First));
501
      end loop;
502
 
503
      return R;
504
   end Vector_Vector_Elementwise_Operation;
505
 
506
   ------------------------------------------------
507
   -- Vector_Vector_Scalar_Elementwise_Operation --
508
   ------------------------------------------------
509
 
510
   function Vector_Vector_Scalar_Elementwise_Operation
511
     (X : X_Vector;
512
      Y : Y_Vector;
513
      Z : Z_Scalar) return Result_Vector
514
   is
515
      R : Result_Vector (X'Range);
516
 
517
   begin
518
      if X'Length /= Y'Length then
519
         raise Constraint_Error with
520
           "vectors are of different length in elementwise operation";
521
      end if;
522
 
523
      for J in R'Range loop
524
         R (J) := Operation (X (J), Y (J - X'First + Y'First), Z);
525
      end loop;
526
 
527
      return R;
528
   end Vector_Vector_Scalar_Elementwise_Operation;
529
 
530
   -----------------------------------------
531
   -- Matrix_Scalar_Elementwise_Operation --
532
   -----------------------------------------
533
 
534
   function Matrix_Scalar_Elementwise_Operation
535
     (Left  : Left_Matrix;
536
      Right : Right_Scalar) return Result_Matrix
537
   is
538
      R : Result_Matrix (Left'Range (1), Left'Range (2));
539
 
540
   begin
541
      for J in R'Range (1) loop
542
         for K in R'Range (2) loop
543
            R (J, K) := Operation (Left (J, K), Right);
544
         end loop;
545
      end loop;
546
 
547
      return R;
548
   end Matrix_Scalar_Elementwise_Operation;
549
 
550
   -----------------------------------------
551
   -- Vector_Scalar_Elementwise_Operation --
552
   -----------------------------------------
553
 
554
   function Vector_Scalar_Elementwise_Operation
555
     (Left  : Left_Vector;
556
      Right : Right_Scalar) return Result_Vector
557
   is
558
      R : Result_Vector (Left'Range);
559
 
560
   begin
561
      for J in R'Range loop
562
         R (J) := Operation (Left (J), Right);
563
      end loop;
564
 
565
      return R;
566
   end Vector_Scalar_Elementwise_Operation;
567
 
568
   -----------------------------------------
569
   -- Scalar_Matrix_Elementwise_Operation --
570
   -----------------------------------------
571
 
572
   function Scalar_Matrix_Elementwise_Operation
573
     (Left  : Left_Scalar;
574
      Right : Right_Matrix) return Result_Matrix
575
   is
576
      R : Result_Matrix (Right'Range (1), Right'Range (2));
577
 
578
   begin
579
      for J in R'Range (1) loop
580
         for K in R'Range (2) loop
581
            R (J, K) := Operation (Left, Right (J, K));
582
         end loop;
583
      end loop;
584
 
585
      return R;
586
   end Scalar_Matrix_Elementwise_Operation;
587
 
588
   -----------------------------------------
589
   -- Scalar_Vector_Elementwise_Operation --
590
   -----------------------------------------
591
 
592
   function Scalar_Vector_Elementwise_Operation
593
     (Left  : Left_Scalar;
594
      Right : Right_Vector) return Result_Vector
595
   is
596
      R : Result_Vector (Right'Range);
597
 
598
   begin
599
      for J in R'Range loop
600
         R (J) := Operation (Left, Right (J));
601
      end loop;
602
 
603
      return R;
604
   end Scalar_Vector_Elementwise_Operation;
605
 
606
   ----------
607
   -- Sqrt --
608
   ----------
609
 
610
   function Sqrt (X : Real'Base) return Real'Base is
611
      Root, Next : Real'Base;
612
 
613
   begin
614
      --  Be defensive: any comparisons with NaN values will yield False.
615
 
616
      if not (X > 0.0) then
617
         if X = 0.0 then
618
            return X;
619
         else
620
            raise Argument_Error;
621
         end if;
622
 
623
      elsif X > Real'Base'Last then
624
 
625
         --  X is infinity, which is its own square root
626
 
627
         return X;
628
      end if;
629
 
630
      --  Compute an initial estimate based on:
631
 
632
      --     X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0),
633
 
634
      --  where M is the mantissa, R is the radix and E the exponent.
635
 
636
      --  By ignoring the mantissa and ignoring the case of an odd
637
      --  exponent, we get a final error that is at most R. In other words,
638
      --  the result has about a single bit precision.
639
 
640
      Root := Real'Base (Real'Machine_Radix) ** (Real'Exponent (X) / 2);
641
 
642
      --  Because of the poor initial estimate, use the Babylonian method of
643
      --  computing the square root, as it is stable for all inputs. Every step
644
      --  will roughly double the precision of the result. Just a few steps
645
      --  suffice in most cases. Eight iterations should give about 2**8 bits
646
      --  of precision.
647
 
648
      for J in 1 .. 8 loop
649
         Next := (Root + X / Root) / 2.0;
650
         exit when Root = Next;
651
         Root := Next;
652
      end loop;
653
 
654
      return Root;
655
   end Sqrt;
656
 
657
   ---------------------------
658
   -- Matrix_Matrix_Product --
659
   ---------------------------
660
 
661
   function Matrix_Matrix_Product
662
     (Left  : Left_Matrix;
663
      Right : Right_Matrix) return Result_Matrix
664
   is
665
      R : Result_Matrix (Left'Range (1), Right'Range (2));
666
 
667
   begin
668
      if Left'Length (2) /= Right'Length (1) then
669
         raise Constraint_Error with
670
           "incompatible dimensions in matrix multiplication";
671
      end if;
672
 
673
      for J in R'Range (1) loop
674
         for K in R'Range (2) loop
675
            declare
676
               S : Result_Scalar := Zero;
677
 
678
            begin
679
               for M in Left'Range (2) loop
680
                  S := S + Left (J, M) *
681
                             Right (M - Left'First (2) + Right'First (1), K);
682
               end loop;
683
 
684
               R (J, K) := S;
685
            end;
686
         end loop;
687
      end loop;
688
 
689
      return R;
690
   end  Matrix_Matrix_Product;
691
 
692
   ----------------------------
693
   -- Matrix_Vector_Solution --
694
   ----------------------------
695
 
696
   function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is
697
      N   : constant Natural := A'Length (1);
698
      MA  : Matrix := A;
699
      MX  : Matrix (A'Range (1), 1 .. 1);
700
      R   : Vector (A'Range (2));
701
      Det : Scalar;
702
 
703
   begin
704
      if A'Length (2) /= N then
705
         raise Constraint_Error with "matrix is not square";
706
      end if;
707
 
708
      if X'Length /= N then
709
         raise Constraint_Error with "incompatible vector length";
710
      end if;
711
 
712
      for J in 0 .. MX'Length (1) - 1 loop
713
         MX (MX'First (1) + J, 1) := X (X'First + J);
714
      end loop;
715
 
716
      Forward_Eliminate (MA, MX, Det);
717
      Back_Substitute (MA, MX);
718
 
719
      for J in 0 .. R'Length - 1 loop
720
         R (R'First + J) := MX (MX'First (1) + J, 1);
721
      end loop;
722
 
723
      return R;
724
   end Matrix_Vector_Solution;
725
 
726
   ----------------------------
727
   -- Matrix_Matrix_Solution --
728
   ----------------------------
729
 
730
   function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is
731
      N   : constant Natural := A'Length (1);
732
      MA  : Matrix (A'Range (2), A'Range (2));
733
      MB  : Matrix (A'Range (2), X'Range (2));
734
      Det : Scalar;
735
 
736
   begin
737
      if A'Length (2) /= N then
738
         raise Constraint_Error with "matrix is not square";
739
      end if;
740
 
741
      if X'Length (1) /= N then
742
         raise Constraint_Error with "matrices have unequal number of rows";
743
      end if;
744
 
745
      for J in 0 .. A'Length (1) - 1 loop
746
         for K in MA'Range (2) loop
747
            MA (MA'First (1) + J, K) := A (A'First (1) + J, K);
748
         end loop;
749
 
750
         for K in MB'Range (2) loop
751
            MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
752
         end loop;
753
      end loop;
754
 
755
      Forward_Eliminate (MA, MB, Det);
756
      Back_Substitute (MA, MB);
757
 
758
      return MB;
759
   end Matrix_Matrix_Solution;
760
 
761
   ---------------------------
762
   -- Matrix_Vector_Product --
763
   ---------------------------
764
 
765
   function Matrix_Vector_Product
766
     (Left  : Matrix;
767
      Right : Right_Vector) return Result_Vector
768
   is
769
      R : Result_Vector (Left'Range (1));
770
 
771
   begin
772
      if Left'Length (2) /= Right'Length then
773
         raise Constraint_Error with
774
            "incompatible dimensions in matrix-vector multiplication";
775
      end if;
776
 
777
      for J in Left'Range (1) loop
778
         declare
779
            S : Result_Scalar := Zero;
780
 
781
         begin
782
            for K in Left'Range (2) loop
783
               S := S + Left (J, K) * Right (K - Left'First (2) + Right'First);
784
            end loop;
785
 
786
            R (J) := S;
787
         end;
788
      end loop;
789
 
790
      return R;
791
   end Matrix_Vector_Product;
792
 
793
   -------------------
794
   -- Outer_Product --
795
   -------------------
796
 
797
   function Outer_Product
798
     (Left  : Left_Vector;
799
      Right : Right_Vector) return Matrix
800
   is
801
      R : Matrix (Left'Range, Right'Range);
802
 
803
   begin
804
      for J in R'Range (1) loop
805
         for K in R'Range (2) loop
806
            R (J, K) := Left (J) * Right (K);
807
         end loop;
808
      end loop;
809
 
810
      return R;
811
   end Outer_Product;
812
 
813
   -----------------
814
   -- Swap_Column --
815
   -----------------
816
 
817
   procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is
818
      Temp : Scalar;
819
   begin
820
      for J in A'Range (1) loop
821
         Temp := A (J, Left);
822
         A (J, Left) := A (J, Right);
823
         A (J, Right) := Temp;
824
      end loop;
825
   end Swap_Column;
826
 
827
   ---------------
828
   -- Transpose --
829
   ---------------
830
 
831
   procedure Transpose (A : Matrix; R : out Matrix) is
832
   begin
833
      for J in R'Range (1) loop
834
         for K in R'Range (2) loop
835
            R (J, K) := A (K - R'First (2) + A'First (1),
836
                           J - R'First (1) + A'First (2));
837
         end loop;
838
      end loop;
839
   end Transpose;
840
 
841
   -------------------------------
842
   -- Update_Matrix_With_Matrix --
843
   -------------------------------
844
 
845
   procedure Update_Matrix_With_Matrix (X : in out X_Matrix; Y : Y_Matrix) is
846
   begin
847
      if X'Length (1) /= Y'Length (1)
848
        or else X'Length (2) /= Y'Length (2)
849
      then
850
         raise Constraint_Error with
851
           "matrices are of different dimension in update operation";
852
      end if;
853
 
854
      for J in X'Range (1) loop
855
         for K in X'Range (2) loop
856
            Update (X (J, K), Y (J - X'First (1) + Y'First (1),
857
                                 K - X'First (2) + Y'First (2)));
858
         end loop;
859
      end loop;
860
   end Update_Matrix_With_Matrix;
861
 
862
   -------------------------------
863
   -- Update_Vector_With_Vector --
864
   -------------------------------
865
 
866
   procedure Update_Vector_With_Vector (X : in out X_Vector; Y : Y_Vector) is
867
   begin
868
      if X'Length /= Y'Length then
869
         raise Constraint_Error with
870
           "vectors are of different length in update operation";
871
      end if;
872
 
873
      for J in X'Range loop
874
         Update (X (J), Y (J - X'First + Y'First));
875
      end loop;
876
   end Update_Vector_With_Vector;
877
 
878
   -----------------
879
   -- Unit_Matrix --
880
   -----------------
881
 
882
   function Unit_Matrix
883
     (Order   : Positive;
884
      First_1 : Integer := 1;
885
      First_2 : Integer := 1) return Matrix
886
   is
887
      R : Matrix (First_1 .. Check_Unit_Last (First_1, Order, First_1),
888
                  First_2 .. Check_Unit_Last (First_2, Order, First_2));
889
 
890
   begin
891
      R := (others => (others => Zero));
892
 
893
      for J in 0 .. Order - 1 loop
894
         R (First_1 + J, First_2 + J) := One;
895
      end loop;
896
 
897
      return R;
898
   end Unit_Matrix;
899
 
900
   -----------------
901
   -- Unit_Vector --
902
   -----------------
903
 
904
   function Unit_Vector
905
     (Index : Integer;
906
      Order : Positive;
907
      First : Integer := 1) return Vector
908
   is
909
      R : Vector (First .. Check_Unit_Last (Index, Order, First));
910
   begin
911
      R := (others => Zero);
912
      R (Index) := One;
913
      return R;
914
   end Unit_Vector;
915
 
916
   ---------------------------
917
   -- Vector_Matrix_Product --
918
   ---------------------------
919
 
920
   function Vector_Matrix_Product
921
     (Left  : Left_Vector;
922
      Right : Matrix) return Result_Vector
923
   is
924
      R : Result_Vector (Right'Range (2));
925
 
926
   begin
927
      if Left'Length /= Right'Length (2) then
928
         raise Constraint_Error with
929
           "incompatible dimensions in vector-matrix multiplication";
930
      end if;
931
 
932
      for J in Right'Range (2) loop
933
         declare
934
            S : Result_Scalar := Zero;
935
 
936
         begin
937
            for K in Right'Range (1) loop
938
               S := S + Left (J - Right'First (1) + Left'First) * Right (K, J);
939
            end loop;
940
 
941
            R (J) := S;
942
         end;
943
      end loop;
944
 
945
      return R;
946
   end Vector_Matrix_Product;
947
 
948
end System.Generic_Array_Operations;

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.