C语言 用非零值初始化双精度数组(BLAS)

7y4bm7vi  于 2023-03-01  发布在  其他
关注(0)|答案(4)|浏览(191)

我分配了一个很大的double向量,比如说有100000个元素。在代码中的某个地方,我想将所有元素设置为一个常量,非零值。如果不对所有元素使用for循环,我怎么能做到这一点呢?我还使用了BLAS包,如果它有用的话。

yvgpqqbh

yvgpqqbh1#

您可以使用std::fill#include <algorithm>):

std::fill(v.begin(), v.end(), 1);

当然,这本质上也只是一个循环。

iyr7buue

iyr7buue2#

从你所说的来看,fill是对的。
请注意,也可以构造一个包含指定值的向量:

std::vector<double> vec(100000, 3.14);

因此,如果“在某个时候”意味着“在施工后立即”,那么就这样做。同样,它意味着你可以这样做:

std::vector<double>(100000, 3.14).swap(vec);

如果“在某个点”表示“在改变大小之后立即”,并且您期望/想要重新分配向量(如果您要使其大于其先前的容量,则使用“期望”;如果您要使其小得多并且想要修剪以保存内存,则使用“想要”),则这可能是有用的。

7gcisfzg

7gcisfzg3#

如果不想循环,则始终使用memset()
memset(myarr, 5, arrsize);,以便用全5填充它。注意隐式转换为unsigned char。
概要

#include <string.h>

 void *
 memset(void *b, int c, size_t len);

描述

The memset() function writes len bytes of value c (converted to an
 unsigned char) to the byte string b.

如果向量很大,你需要它快速运行,并且你使用gcc,那么:
块移动(memcpy)和块设置(memset)的代码生成被重写。GCC现在可以根据要复制的块的大小和要优化的CPU选择最佳算法(循环、展开循环、带rep前缀的指令或库调用)。

cczfrluj

cczfrluj4#

不幸的是,其他答案没有按照提示进行,OP希望将数组元素设置为零。使用BLAS而不是更惯用的函数(如memsetfill)可能有几个原因。例如,BLAS操作可以线程化。此外,memsetfill不提供开箱即用的 strided 操作。
乍一看,BLAS库似乎不提供这样的功能,但是有两个选项:
1.为此,可以利用BLAS函数xSCAL(对于不同的数值类型,x可以是sdcz)。
SCAL执行缩放操作V <- a*V。对于a = 0,它将元素设置为零(大部分)。
1.使用xCOPY,并从堆栈内存中复制一个零。
完整代码如下。
这两种方法都有各自的问题,第一种策略依赖于任意浮点数x的任意x*0.0 == 0.0。从技术上讲,这对于x == NANx == infinity是不正确的(两种情况都示出)。也许BLAS可以以实际上给出这一点的非IEEE兼容方式来编译。在任何情况下,如果由于某种原因你知道原始值是正则数,那么你可以使用这个。2另一个问题是你可能得到有符号的零(-0.0)作为元素,这是可以的,除了它们可能最终不都是相同的均匀(例如,正零,0.0)。
第二种更健壮,但是依赖于BLAS接受零步长值。(BLAS是在70年代用Fortran发明和编码的,当时整数零还没有发明。)我知道的大多数实现都允许零增量,至少对于xCOPY是这样。它还需要从某个地方取出“第一个”零;在本例中,仅在堆栈中创建。(如果要泛化到GPU BLAS(cuBLAS),则需要在GPU中分配此零。)
因此,换句话说,您必须了解您的平台和可用的BLAS。

#include<cstdint>
#include<iostream>
#include<limits>

extern "C" {
void sscal_(int32_t const& n, float const& a, float* x, int32_t const& incx);
void scopy_(int32_t const& n, float const* x, int32_t const& incx, float* y, int32_t const& incy);
}

void set_zero_1(int32_t n, float* x, int incx) {
  sscal_(n, 0.0F, x, incx);
}

void fill_value(int32_t n, float* x, int incx, float const* value_ptr) {
  scopy_(n, value_ptr, 0, x, incx);
}

void set_zero_2(int32_t n, float* x, int incx) {
  float const value = 0.0F;  // can also be allocated or be a global if necessary
  fill_value(n, x, incx, &value);
}

int main() {
  float X[12] = {
    99.9, 0.0, 0.0, 
    std::numeric_limits<float>::quiet_NaN(), 0.0, 0.0, 
    std::numeric_limits<float>::infinity(), 0.0, 0.0, 
    99.9, 0.0, 0.0
  };
//set_zero_1(                 4,            &X[0],            3);  // this fails because NAN and INF
  set_zero_2(/*num elements*/ 4, /*origin*/ &X[0], /*stride*/ 3);
  for(int i = 0; i != 12; i += 3) std::cout << X[i] << std::endl;  // prints zeros
}

像这样使用。它将打印零。

$ c++ a.cpp -L/opt/intel/oneapi/mkl/2023.0.0/lib/intel64 -lmkl_rt 
$ export LD_LIBRARY_PATH=/opt/intel/oneapi/mkl/2023.0.0/lib/intel64
$ ./a.out
0
0
0
0

另外,也有这种方法,但它依赖于BLAS是顺序的,而且它不能并行化或矢量化任何操作,所以这可能是最糟糕的选择。

void set_zero_3(int32_t n, float* x, int incx) {
  *x = 0.0F;  // set one element to zero, somehow. (see note above about memory that is not accessible from the CPU)
  scopy_(n - 1, x, incx, x + incx, incx);
}

相关问题