1 //===-- runtime/transformational.cpp --------------------------------------===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 
9 // Implements the transformational intrinsic functions of Fortran 2018 that
10 // rearrange or duplicate data without (much) regard to type.  These are
11 // CSHIFT, EOSHIFT, PACK, RESHAPE, SPREAD, TRANSPOSE, and UNPACK.
12 //
13 // Many of these are defined in the 2018 standard with text that makes sense
14 // only if argument arrays have lower bounds of one.  Rather than interpret
15 // these cases as implying a hidden constraint, these implementations
16 // work with arbitrary lower bounds.  This may be technically an extension
17 // of the standard but it more likely to conform with its intent.
18 
19 #include "flang/Runtime/transformational.h"
20 #include "copy.h"
21 #include "terminator.h"
22 #include "tools.h"
23 #include "flang/Runtime/descriptor.h"
24 #include <algorithm>
25 
26 namespace Fortran::runtime {
27 
28 // Utility for CSHIFT & EOSHIFT rank > 1 cases that determines the shift count
29 // for each of the vector sections of the result.
30 class ShiftControl {
31 public:
ShiftControl(const Descriptor & s,Terminator & t,int dim)32   ShiftControl(const Descriptor &s, Terminator &t, int dim)
33       : shift_{s}, terminator_{t}, shiftRank_{s.rank()}, dim_{dim} {}
Init(const Descriptor & source,const char * which)34   void Init(const Descriptor &source, const char *which) {
35     int rank{source.rank()};
36     RUNTIME_CHECK(terminator_, shiftRank_ == 0 || shiftRank_ == rank - 1);
37     auto catAndKind{shift_.type().GetCategoryAndKind()};
38     RUNTIME_CHECK(
39         terminator_, catAndKind && catAndKind->first == TypeCategory::Integer);
40     shiftElemLen_ = catAndKind->second;
41     if (shiftRank_ > 0) {
42       int k{0};
43       for (int j{0}; j < rank; ++j) {
44         if (j + 1 != dim_) {
45           const Dimension &shiftDim{shift_.GetDimension(k)};
46           lb_[k++] = shiftDim.LowerBound();
47           if (shiftDim.Extent() != source.GetDimension(j).Extent()) {
48             terminator_.Crash("%s: on dimension %d, SHIFT= has extent %jd but "
49                               "SOURCE= has extent %jd",
50                 which, k, static_cast<std::intmax_t>(shiftDim.Extent()),
51                 static_cast<std::intmax_t>(source.GetDimension(j).Extent()));
52           }
53         }
54       }
55     } else {
56       shiftCount_ =
57           GetInt64(shift_.OffsetElement<char>(), shiftElemLen_, terminator_);
58     }
59   }
GetShift(const SubscriptValue resultAt[]) const60   SubscriptValue GetShift(const SubscriptValue resultAt[]) const {
61     if (shiftRank_ > 0) {
62       SubscriptValue shiftAt[maxRank];
63       int k{0};
64       for (int j{0}; j < shiftRank_ + 1; ++j) {
65         if (j + 1 != dim_) {
66           shiftAt[k] = lb_[k] + resultAt[j] - 1;
67           ++k;
68         }
69       }
70       return GetInt64(
71           shift_.Element<char>(shiftAt), shiftElemLen_, terminator_);
72     } else {
73       return shiftCount_; // invariant count extracted in Init()
74     }
75   }
76 
77 private:
78   const Descriptor &shift_;
79   Terminator &terminator_;
80   int shiftRank_;
81   int dim_;
82   SubscriptValue lb_[maxRank];
83   std::size_t shiftElemLen_;
84   SubscriptValue shiftCount_{};
85 };
86 
87 // Fill an EOSHIFT result with default boundary values
DefaultInitialize(const Descriptor & result,Terminator & terminator)88 static void DefaultInitialize(
89     const Descriptor &result, Terminator &terminator) {
90   auto catAndKind{result.type().GetCategoryAndKind()};
91   RUNTIME_CHECK(
92       terminator, catAndKind && catAndKind->first != TypeCategory::Derived);
93   std::size_t elementLen{result.ElementBytes()};
94   std::size_t bytes{result.Elements() * elementLen};
95   if (catAndKind->first == TypeCategory::Character) {
96     switch (int kind{catAndKind->second}) {
97     case 1:
98       std::fill_n(result.OffsetElement<char>(), bytes, ' ');
99       break;
100     case 2:
101       std::fill_n(result.OffsetElement<char16_t>(), bytes / 2,
102           static_cast<char16_t>(' '));
103       break;
104     case 4:
105       std::fill_n(result.OffsetElement<char32_t>(), bytes / 4,
106           static_cast<char32_t>(' '));
107       break;
108     default:
109       terminator.Crash("not yet implemented: EOSHIFT: CHARACTER kind %d", kind);
110     }
111   } else {
112     std::memset(result.raw().base_addr, 0, bytes);
113   }
114 }
115 
AllocateResult(Descriptor & result,const Descriptor & source,int rank,const SubscriptValue extent[],Terminator & terminator,const char * function)116 static inline std::size_t AllocateResult(Descriptor &result,
117     const Descriptor &source, int rank, const SubscriptValue extent[],
118     Terminator &terminator, const char *function) {
119   std::size_t elementLen{source.ElementBytes()};
120   const DescriptorAddendum *sourceAddendum{source.Addendum()};
121   result.Establish(source.type(), elementLen, nullptr, rank, extent,
122       CFI_attribute_allocatable, sourceAddendum != nullptr);
123   if (sourceAddendum) {
124     *result.Addendum() = *sourceAddendum;
125   }
126   for (int j{0}; j < rank; ++j) {
127     result.GetDimension(j).SetBounds(1, extent[j]);
128   }
129   if (int stat{result.Allocate()}) {
130     terminator.Crash(
131         "%s: Could not allocate memory for result (stat=%d)", function, stat);
132   }
133   return elementLen;
134 }
135 
136 extern "C" {
137 
138 // CSHIFT where rank of ARRAY argument > 1
RTNAME(Cshift)139 void RTNAME(Cshift)(Descriptor &result, const Descriptor &source,
140     const Descriptor &shift, int dim, const char *sourceFile, int line) {
141   Terminator terminator{sourceFile, line};
142   int rank{source.rank()};
143   RUNTIME_CHECK(terminator, rank > 1);
144   if (dim < 1 || dim > rank) {
145     terminator.Crash(
146         "CSHIFT: DIM=%d must be >= 1 and <= SOURCE= rank %d", dim, rank);
147   }
148   ShiftControl shiftControl{shift, terminator, dim};
149   shiftControl.Init(source, "CSHIFT");
150   SubscriptValue extent[maxRank];
151   source.GetShape(extent);
152   AllocateResult(result, source, rank, extent, terminator, "CSHIFT");
153   SubscriptValue resultAt[maxRank];
154   for (int j{0}; j < rank; ++j) {
155     resultAt[j] = 1;
156   }
157   SubscriptValue sourceLB[maxRank];
158   source.GetLowerBounds(sourceLB);
159   SubscriptValue dimExtent{extent[dim - 1]};
160   SubscriptValue dimLB{sourceLB[dim - 1]};
161   SubscriptValue &resDim{resultAt[dim - 1]};
162   for (std::size_t n{result.Elements()}; n > 0; n -= dimExtent) {
163     SubscriptValue shiftCount{shiftControl.GetShift(resultAt)};
164     SubscriptValue sourceAt[maxRank];
165     for (int j{0}; j < rank; ++j) {
166       sourceAt[j] = sourceLB[j] + resultAt[j] - 1;
167     }
168     SubscriptValue &sourceDim{sourceAt[dim - 1]};
169     sourceDim = dimLB + shiftCount % dimExtent;
170     if (sourceDim < dimLB) {
171       sourceDim += dimExtent;
172     }
173     for (resDim = 1; resDim <= dimExtent; ++resDim) {
174       CopyElement(result, resultAt, source, sourceAt, terminator);
175       if (++sourceDim == dimLB + dimExtent) {
176         sourceDim = dimLB;
177       }
178     }
179     result.IncrementSubscripts(resultAt);
180   }
181 }
182 
183 // CSHIFT where rank of ARRAY argument == 1
RTNAME(CshiftVector)184 void RTNAME(CshiftVector)(Descriptor &result, const Descriptor &source,
185     std::int64_t shift, const char *sourceFile, int line) {
186   Terminator terminator{sourceFile, line};
187   RUNTIME_CHECK(terminator, source.rank() == 1);
188   const Dimension &sourceDim{source.GetDimension(0)};
189   SubscriptValue extent{sourceDim.Extent()};
190   AllocateResult(result, source, 1, &extent, terminator, "CSHIFT");
191   SubscriptValue lb{sourceDim.LowerBound()};
192   for (SubscriptValue j{0}; j < extent; ++j) {
193     SubscriptValue resultAt{1 + j};
194     SubscriptValue sourceAt{lb + (j + shift) % extent};
195     if (sourceAt < lb) {
196       sourceAt += extent;
197     }
198     CopyElement(result, &resultAt, source, &sourceAt, terminator);
199   }
200 }
201 
202 // EOSHIFT of rank > 1
RTNAME(Eoshift)203 void RTNAME(Eoshift)(Descriptor &result, const Descriptor &source,
204     const Descriptor &shift, const Descriptor *boundary, int dim,
205     const char *sourceFile, int line) {
206   Terminator terminator{sourceFile, line};
207   SubscriptValue extent[maxRank];
208   int rank{source.GetShape(extent)};
209   RUNTIME_CHECK(terminator, rank > 1);
210   if (dim < 1 || dim > rank) {
211     terminator.Crash(
212         "EOSHIFT: DIM=%d must be >= 1 and <= SOURCE= rank %d", dim, rank);
213   }
214   std::size_t elementLen{
215       AllocateResult(result, source, rank, extent, terminator, "EOSHIFT")};
216   int boundaryRank{-1};
217   if (boundary) {
218     boundaryRank = boundary->rank();
219     RUNTIME_CHECK(terminator, boundaryRank == 0 || boundaryRank == rank - 1);
220     RUNTIME_CHECK(terminator, boundary->type() == source.type());
221     if (boundary->ElementBytes() != elementLen) {
222       terminator.Crash("EOSHIFT: BOUNDARY= has element byte length %zd, but "
223                        "SOURCE= has length %zd",
224           boundary->ElementBytes(), elementLen);
225     }
226     if (boundaryRank > 0) {
227       int k{0};
228       for (int j{0}; j < rank; ++j) {
229         if (j != dim - 1) {
230           if (boundary->GetDimension(k).Extent() != extent[j]) {
231             terminator.Crash("EOSHIFT: BOUNDARY= has extent %jd on dimension "
232                              "%d but must conform with extent %jd of SOURCE=",
233                 static_cast<std::intmax_t>(boundary->GetDimension(k).Extent()),
234                 k + 1, static_cast<std::intmax_t>(extent[j]));
235           }
236           ++k;
237         }
238       }
239     }
240   }
241   ShiftControl shiftControl{shift, terminator, dim};
242   shiftControl.Init(source, "EOSHIFT");
243   SubscriptValue resultAt[maxRank];
244   for (int j{0}; j < rank; ++j) {
245     resultAt[j] = 1;
246   }
247   if (!boundary) {
248     DefaultInitialize(result, terminator);
249   }
250   SubscriptValue sourceLB[maxRank];
251   source.GetLowerBounds(sourceLB);
252   SubscriptValue boundaryAt[maxRank];
253   if (boundaryRank > 0) {
254     boundary->GetLowerBounds(boundaryAt);
255   }
256   SubscriptValue dimExtent{extent[dim - 1]};
257   SubscriptValue dimLB{sourceLB[dim - 1]};
258   SubscriptValue &resDim{resultAt[dim - 1]};
259   for (std::size_t n{result.Elements()}; n > 0; n -= dimExtent) {
260     SubscriptValue shiftCount{shiftControl.GetShift(resultAt)};
261     SubscriptValue sourceAt[maxRank];
262     for (int j{0}; j < rank; ++j) {
263       sourceAt[j] = sourceLB[j] + resultAt[j] - 1;
264     }
265     SubscriptValue &sourceDim{sourceAt[dim - 1]};
266     sourceDim = dimLB + shiftCount;
267     for (resDim = 1; resDim <= dimExtent; ++resDim) {
268       if (sourceDim >= dimLB && sourceDim < dimLB + dimExtent) {
269         CopyElement(result, resultAt, source, sourceAt, terminator);
270       } else if (boundary) {
271         CopyElement(result, resultAt, *boundary, boundaryAt, terminator);
272       }
273       ++sourceDim;
274     }
275     result.IncrementSubscripts(resultAt);
276     if (boundaryRank > 0) {
277       boundary->IncrementSubscripts(boundaryAt);
278     }
279   }
280 }
281 
282 // EOSHIFT of vector
RTNAME(EoshiftVector)283 void RTNAME(EoshiftVector)(Descriptor &result, const Descriptor &source,
284     std::int64_t shift, const Descriptor *boundary, const char *sourceFile,
285     int line) {
286   Terminator terminator{sourceFile, line};
287   RUNTIME_CHECK(terminator, source.rank() == 1);
288   SubscriptValue extent{source.GetDimension(0).Extent()};
289   std::size_t elementLen{
290       AllocateResult(result, source, 1, &extent, terminator, "EOSHIFT")};
291   if (boundary) {
292     RUNTIME_CHECK(terminator, boundary->rank() == 0);
293     RUNTIME_CHECK(terminator, boundary->type() == source.type());
294     if (boundary->ElementBytes() != elementLen) {
295       terminator.Crash("EOSHIFT: BOUNDARY= has element byte length %zd but "
296                        "SOURCE= has length %zd",
297           boundary->ElementBytes(), elementLen);
298     }
299   }
300   if (!boundary) {
301     DefaultInitialize(result, terminator);
302   }
303   SubscriptValue lb{source.GetDimension(0).LowerBound()};
304   for (SubscriptValue j{1}; j <= extent; ++j) {
305     SubscriptValue sourceAt{lb + j - 1 + shift};
306     if (sourceAt >= lb && sourceAt < lb + extent) {
307       CopyElement(result, &j, source, &sourceAt, terminator);
308     } else if (boundary) {
309       CopyElement(result, &j, *boundary, 0, terminator);
310     }
311   }
312 }
313 
314 // PACK
RTNAME(Pack)315 void RTNAME(Pack)(Descriptor &result, const Descriptor &source,
316     const Descriptor &mask, const Descriptor *vector, const char *sourceFile,
317     int line) {
318   Terminator terminator{sourceFile, line};
319   CheckConformability(source, mask, terminator, "PACK", "ARRAY=", "MASK=");
320   auto maskType{mask.type().GetCategoryAndKind()};
321   RUNTIME_CHECK(
322       terminator, maskType && maskType->first == TypeCategory::Logical);
323   SubscriptValue trues{0};
324   if (mask.rank() == 0) {
325     if (IsLogicalElementTrue(mask, nullptr)) {
326       trues = source.Elements();
327     }
328   } else {
329     SubscriptValue maskAt[maxRank];
330     mask.GetLowerBounds(maskAt);
331     for (std::size_t n{mask.Elements()}; n > 0; --n) {
332       if (IsLogicalElementTrue(mask, maskAt)) {
333         ++trues;
334       }
335       mask.IncrementSubscripts(maskAt);
336     }
337   }
338   SubscriptValue extent{trues};
339   if (vector) {
340     RUNTIME_CHECK(terminator, vector->rank() == 1);
341     RUNTIME_CHECK(terminator, source.type() == vector->type());
342     if (source.ElementBytes() != vector->ElementBytes()) {
343       terminator.Crash("PACK: SOURCE= has element byte length %zd, but VECTOR= "
344                        "has length %zd",
345           source.ElementBytes(), vector->ElementBytes());
346     }
347     extent = vector->GetDimension(0).Extent();
348     if (extent < trues) {
349       terminator.Crash("PACK: VECTOR= has extent %jd but there are %jd MASK= "
350                        "elements that are .TRUE.",
351           static_cast<std::intmax_t>(extent),
352           static_cast<std::intmax_t>(trues));
353     }
354   }
355   AllocateResult(result, source, 1, &extent, terminator, "PACK");
356   SubscriptValue sourceAt[maxRank], resultAt{1};
357   source.GetLowerBounds(sourceAt);
358   if (mask.rank() == 0) {
359     if (IsLogicalElementTrue(mask, nullptr)) {
360       for (SubscriptValue n{trues}; n > 0; --n) {
361         CopyElement(result, &resultAt, source, sourceAt, terminator);
362         ++resultAt;
363         source.IncrementSubscripts(sourceAt);
364       }
365     }
366   } else {
367     SubscriptValue maskAt[maxRank];
368     mask.GetLowerBounds(maskAt);
369     for (std::size_t n{source.Elements()}; n > 0; --n) {
370       if (IsLogicalElementTrue(mask, maskAt)) {
371         CopyElement(result, &resultAt, source, sourceAt, terminator);
372         ++resultAt;
373       }
374       source.IncrementSubscripts(sourceAt);
375       mask.IncrementSubscripts(maskAt);
376     }
377   }
378   if (vector) {
379     SubscriptValue vectorAt{
380         vector->GetDimension(0).LowerBound() + resultAt - 1};
381     for (; resultAt <= extent; ++resultAt, ++vectorAt) {
382       CopyElement(result, &resultAt, *vector, &vectorAt, terminator);
383     }
384   }
385 }
386 
387 // RESHAPE
388 // F2018 16.9.163
RTNAME(Reshape)389 void RTNAME(Reshape)(Descriptor &result, const Descriptor &source,
390     const Descriptor &shape, const Descriptor *pad, const Descriptor *order,
391     const char *sourceFile, int line) {
392   // Compute and check the rank of the result.
393   Terminator terminator{sourceFile, line};
394   RUNTIME_CHECK(terminator, shape.rank() == 1);
395   RUNTIME_CHECK(terminator, shape.type().IsInteger());
396   SubscriptValue resultRank{shape.GetDimension(0).Extent()};
397   if (resultRank < 0 || resultRank > static_cast<SubscriptValue>(maxRank)) {
398     terminator.Crash(
399         "RESHAPE: SHAPE= vector length %jd implies a bad result rank",
400         static_cast<std::intmax_t>(resultRank));
401   }
402 
403   // Extract and check the shape of the result; compute its element count.
404   SubscriptValue resultExtent[maxRank];
405   std::size_t shapeElementBytes{shape.ElementBytes()};
406   std::size_t resultElements{1};
407   SubscriptValue shapeSubscript{shape.GetDimension(0).LowerBound()};
408   for (int j{0}; j < resultRank; ++j, ++shapeSubscript) {
409     resultExtent[j] = GetInt64(
410         shape.Element<char>(&shapeSubscript), shapeElementBytes, terminator);
411     if (resultExtent[j] < 0) {
412       terminator.Crash("RESHAPE: bad value for SHAPE(%d)=%jd", j + 1,
413           static_cast<std::intmax_t>(resultExtent[j]));
414     }
415     resultElements *= resultExtent[j];
416   }
417 
418   // Check that there are sufficient elements in the SOURCE=, or that
419   // the optional PAD= argument is present and nonempty.
420   std::size_t elementBytes{source.ElementBytes()};
421   std::size_t sourceElements{source.Elements()};
422   std::size_t padElements{pad ? pad->Elements() : 0};
423   if (resultElements > sourceElements) {
424     if (padElements <= 0) {
425       terminator.Crash(
426           "RESHAPE: not enough elements, need %zd but only have %zd",
427           resultElements, sourceElements);
428     }
429     if (pad->ElementBytes() != elementBytes) {
430       terminator.Crash("RESHAPE: PAD= has element byte length %zd but SOURCE= "
431                        "has length %zd",
432           pad->ElementBytes(), elementBytes);
433     }
434   }
435 
436   // Extract and check the optional ORDER= argument, which must be a
437   // permutation of [1..resultRank].
438   int dimOrder[maxRank];
439   if (order) {
440     RUNTIME_CHECK(terminator, order->rank() == 1);
441     RUNTIME_CHECK(terminator, order->type().IsInteger());
442     if (order->GetDimension(0).Extent() != resultRank) {
443       terminator.Crash("RESHAPE: the extent of ORDER (%jd) must match the rank"
444                        " of the SHAPE (%d)",
445           static_cast<std::intmax_t>(order->GetDimension(0).Extent()),
446           resultRank);
447     }
448     std::uint64_t values{0};
449     SubscriptValue orderSubscript{order->GetDimension(0).LowerBound()};
450     std::size_t orderElementBytes{order->ElementBytes()};
451     for (SubscriptValue j{0}; j < resultRank; ++j, ++orderSubscript) {
452       auto k{GetInt64(order->Element<char>(&orderSubscript), orderElementBytes,
453           terminator)};
454       if (k < 1 || k > resultRank || ((values >> k) & 1)) {
455         terminator.Crash("RESHAPE: bad value for ORDER element (%jd)",
456             static_cast<std::intmax_t>(k));
457       }
458       values |= std::uint64_t{1} << k;
459       dimOrder[j] = k - 1;
460     }
461   } else {
462     for (int j{0}; j < resultRank; ++j) {
463       dimOrder[j] = j;
464     }
465   }
466 
467   // Allocate result descriptor
468   AllocateResult(
469       result, source, resultRank, resultExtent, terminator, "RESHAPE");
470 
471   // Populate the result's elements.
472   SubscriptValue resultSubscript[maxRank];
473   result.GetLowerBounds(resultSubscript);
474   SubscriptValue sourceSubscript[maxRank];
475   source.GetLowerBounds(sourceSubscript);
476   std::size_t resultElement{0};
477   std::size_t elementsFromSource{std::min(resultElements, sourceElements)};
478   for (; resultElement < elementsFromSource; ++resultElement) {
479     CopyElement(result, resultSubscript, source, sourceSubscript, terminator);
480     source.IncrementSubscripts(sourceSubscript);
481     result.IncrementSubscripts(resultSubscript, dimOrder);
482   }
483   if (resultElement < resultElements) {
484     // Remaining elements come from the optional PAD= argument.
485     SubscriptValue padSubscript[maxRank];
486     pad->GetLowerBounds(padSubscript);
487     for (; resultElement < resultElements; ++resultElement) {
488       CopyElement(result, resultSubscript, *pad, padSubscript, terminator);
489       pad->IncrementSubscripts(padSubscript);
490       result.IncrementSubscripts(resultSubscript, dimOrder);
491     }
492   }
493 }
494 
495 // SPREAD
RTNAME(Spread)496 void RTNAME(Spread)(Descriptor &result, const Descriptor &source, int dim,
497     std::int64_t ncopies, const char *sourceFile, int line) {
498   Terminator terminator{sourceFile, line};
499   int rank{source.rank() + 1};
500   RUNTIME_CHECK(terminator, rank <= maxRank);
501   if (dim < 1 || dim > rank) {
502     terminator.Crash("SPREAD: DIM=%d argument for rank-%d source array "
503                      "must be greater than 1 and less than or equal to %d",
504         dim, rank - 1, rank);
505   }
506   ncopies = std::max<std::int64_t>(ncopies, 0);
507   SubscriptValue extent[maxRank];
508   int k{0};
509   for (int j{0}; j < rank; ++j) {
510     extent[j] = j == dim - 1 ? ncopies : source.GetDimension(k++).Extent();
511   }
512   AllocateResult(result, source, rank, extent, terminator, "SPREAD");
513   SubscriptValue resultAt[maxRank];
514   for (int j{0}; j < rank; ++j) {
515     resultAt[j] = 1;
516   }
517   SubscriptValue &resultDim{resultAt[dim - 1]};
518   SubscriptValue sourceAt[maxRank];
519   source.GetLowerBounds(sourceAt);
520   for (std::size_t n{result.Elements()}; n > 0; n -= ncopies) {
521     for (resultDim = 1; resultDim <= ncopies; ++resultDim) {
522       CopyElement(result, resultAt, source, sourceAt, terminator);
523     }
524     result.IncrementSubscripts(resultAt);
525     source.IncrementSubscripts(sourceAt);
526   }
527 }
528 
529 // TRANSPOSE
RTNAME(Transpose)530 void RTNAME(Transpose)(Descriptor &result, const Descriptor &matrix,
531     const char *sourceFile, int line) {
532   Terminator terminator{sourceFile, line};
533   RUNTIME_CHECK(terminator, matrix.rank() == 2);
534   SubscriptValue extent[2]{
535       matrix.GetDimension(1).Extent(), matrix.GetDimension(0).Extent()};
536   AllocateResult(result, matrix, 2, extent, terminator, "TRANSPOSE");
537   SubscriptValue resultAt[2]{1, 1};
538   SubscriptValue matrixLB[2];
539   matrix.GetLowerBounds(matrixLB);
540   for (std::size_t n{result.Elements()}; n-- > 0;
541        result.IncrementSubscripts(resultAt)) {
542     SubscriptValue matrixAt[2]{
543         matrixLB[0] + resultAt[1] - 1, matrixLB[1] + resultAt[0] - 1};
544     CopyElement(result, resultAt, matrix, matrixAt, terminator);
545   }
546 }
547 
548 // UNPACK
RTNAME(Unpack)549 void RTNAME(Unpack)(Descriptor &result, const Descriptor &vector,
550     const Descriptor &mask, const Descriptor &field, const char *sourceFile,
551     int line) {
552   Terminator terminator{sourceFile, line};
553   RUNTIME_CHECK(terminator, vector.rank() == 1);
554   int rank{mask.rank()};
555   RUNTIME_CHECK(terminator, rank > 0);
556   SubscriptValue extent[maxRank];
557   mask.GetShape(extent);
558   CheckConformability(mask, field, terminator, "UNPACK", "MASK=", "FIELD=");
559   std::size_t elementLen{
560       AllocateResult(result, field, rank, extent, terminator, "UNPACK")};
561   RUNTIME_CHECK(terminator, vector.type() == field.type());
562   if (vector.ElementBytes() != elementLen) {
563     terminator.Crash(
564         "UNPACK: VECTOR= has element byte length %zd but FIELD= has length %zd",
565         vector.ElementBytes(), elementLen);
566   }
567   SubscriptValue resultAt[maxRank], maskAt[maxRank], fieldAt[maxRank],
568       vectorAt{vector.GetDimension(0).LowerBound()};
569   for (int j{0}; j < rank; ++j) {
570     resultAt[j] = 1;
571   }
572   mask.GetLowerBounds(maskAt);
573   field.GetLowerBounds(fieldAt);
574   SubscriptValue vectorElements{vector.GetDimension(0).Extent()};
575   SubscriptValue vectorLeft{vectorElements};
576   for (std::size_t n{result.Elements()}; n-- > 0;) {
577     if (IsLogicalElementTrue(mask, maskAt)) {
578       if (vectorLeft-- == 0) {
579         terminator.Crash(
580             "UNPACK: VECTOR= argument has fewer elements (%d) than "
581             "MASK= has .TRUE. entries",
582             vectorElements);
583       }
584       CopyElement(result, resultAt, vector, &vectorAt, terminator);
585       ++vectorAt;
586     } else {
587       CopyElement(result, resultAt, field, fieldAt, terminator);
588     }
589     result.IncrementSubscripts(resultAt);
590     mask.IncrementSubscripts(maskAt);
591     field.IncrementSubscripts(fieldAt);
592   }
593 }
594 
595 } // extern "C"
596 } // namespace Fortran::runtime
597