NumPy 源码解析（二十六）
.\numpy\numpy\linalg\lapack_lite\fortran.py
# WARNING! This a Python 2 script. Read README.rst for rationale.
# 引入 re 和 itertools 模块
import re
import itertools
# 检查给定行是否为空白行
def isBlank(line):
return not line
# 检查给定行是否是标签行（以数字开头）
def isLabel(line):
return line[0].isdigit()
# 检查给定行是否是注释行（不以空格开头）
def isComment(line):
return line[0] != ' '
# 检查给定行是否是续行（第 6 个字符不是空格）
def isContinuation(line):
return line[5] != ' '
# 定义常量 COMMENT, STATEMENT, CONTINUATION 用于表示不同类型的行
COMMENT, STATEMENT, CONTINUATION = 0, 1, 2
# 函数：确定 Fortran 代码行的类型
def lineType(line):
"""Return the type of a line of Fortran code."""
if isBlank(line):
return COMMENT
elif isLabel(line):
return STATEMENT
elif isComment(line):
return COMMENT
elif isContinuation(line):
return CONTINUATION
else:
return STATEMENT
# 类：LineIterator，用于迭代处理行并去除行尾空格
class LineIterator:
"""LineIterator(iterable)
Return rstrip()'d lines from iterable, while keeping a count of the
line number in the .lineno attribute.
"""
def __init__(self, iterable):
object.__init__(self)
self.iterable = iter(iterable)
self.lineno = 0
def __iter__(self):
return self
def __next__(self):
self.lineno += 1
line = next(self.iterable)
line = line.rstrip()
return line
next = __next__
# 类：PushbackIterator，支持将元素推回迭代器
class PushbackIterator:
"""PushbackIterator(iterable)
Return an iterator for which items can be pushed back into.
Call the .pushback(item) method to have item returned as the next
value of next().
"""
def __init__(self, iterable):
object.__init__(self)
self.iterable = iter(iterable)
self.buffer = []
def __iter__(self):
return self
def __next__(self):
if self.buffer:
return self.buffer.pop()
else:
return next(self.iterable)
def pushback(self, item):
self.buffer.append(item)
next = __next__
# 函数：fortranSourceLines，返回一个迭代器，处理 Fortran 源文件的语句行
def fortranSourceLines(fo):
"""Return an iterator over statement lines of a Fortran source file.
Comment and blank lines are stripped out, and continuation lines are
merged.
"""
numberingiter = LineIterator(fo)
# 在末尾添加一个额外的空字符串，以处理续行
with_extra = itertools.chain(numberingiter, [''])
pushbackiter = PushbackIterator(with_extra)
for line in pushbackiter:
t = lineType(line)
if t == COMMENT:
continue
elif t == STATEMENT:
lines = [line]
# 处理续行的逻辑，保证行的连续性
for next_line in pushbackiter:
t = lineType(next_line)
if t == CONTINUATION:
lines.append(next_line[6:])
else:
pushbackiter.pushback(next_line)
break
yield numberingiter.lineno, ''.join(lines)
else:
raise ValueError("jammed: continuation line not expected: %s:%d" %
(fo.name, numberingiter.lineno))
# 函数：getDependencies(filename)
def getDependencies(filename):
"""
对于一个 Fortran 源文件，返回其中声明为 EXTERNAL 的例程列表。
"""
# 编译正则表达式模式，用于匹配以 EXTERNAL 开头的行（忽略大小写）
external_pat = re.compile(r'^\s*EXTERNAL\s', re.I)
# 初始化例程列表
routines = []
# 打开文件并进行迭代处理每一行
with open(filename) as fo:
# 使用自定义的函数迭代处理 Fortran 源文件的行
for lineno, line in fortranSourceLines(fo):
# 尝试在当前行中匹配 EXTERNAL 的模式
m = external_pat.match(line)
if m:
# 如果匹配成功，提取 EXTERNAL 后面的例程名称列表
names = line[m.end():].strip().split(',')
# 去除每个名称的首尾空白字符并转换为小写
names = [n.strip().lower() for n in names]
# 过滤掉空的名称
names = [n for n in names if n]
# 将提取的例程名称列表添加到总例程列表中
routines.extend(names)
# 返回最终的例程名称列表作为结果
return routines
.\numpy\numpy\linalg\lapack_lite\lapack_lite_names.h
/*
* NOTE: This is generated code. Look in numpy/linalg/lapack_lite for
* information on remaking this file.
*/
/*
* This file renames all BLAS/LAPACK and f2c symbols to avoid
* dynamic symbol name conflicts, in cases where e.g.
* integer sizes do not match with 'standard' ABI.
*/
// 定义了一系列宏，用于重命名 BLAS/LAPACK 和 f2c 符号，避免动态符号名称冲突
#define caxpy_ BLAS_FUNC(caxpy)
#define ccopy_ BLAS_FUNC(ccopy)
#define cdotc_ BLAS_FUNC(cdotc)
#define cdotu_ BLAS_FUNC(cdotu)
#define cgebak_ BLAS_FUNC(cgebak)
#define cgebal_ BLAS_FUNC(cgebal)
#define cgebd2_ BLAS_FUNC(cgebd2)
#define cgebrd_ BLAS_FUNC(cgebrd)
#define cgeev_ BLAS_FUNC(cgeev)
#define cgehd2_ BLAS_FUNC(cgehd2)
#define cgehrd_ BLAS_FUNC(cgehrd)
#define cgelq2_ BLAS_FUNC(cgelq2)
#define cgelqf_ BLAS_FUNC(cgelqf)
#define cgelsd_ BLAS_FUNC(cgelsd)
#define cgemm_ BLAS_FUNC(cgemm)
#define cgemv_ BLAS_FUNC(cgemv)
#define cgeqr2_ BLAS_FUNC(cgeqr2)
#define cgeqrf_ BLAS_FUNC(cgeqrf)
#define cgerc_ BLAS_FUNC(cgerc)
#define cgeru_ BLAS_FUNC(cgeru)
#define cgesdd_ BLAS_FUNC(cgesdd)
#define cgesv_ BLAS_FUNC(cgesv)
#define cgetf2_ BLAS_FUNC(cgetf2)
#define cgetrf_ BLAS_FUNC(cgetrf)
#define cgetrs_ BLAS_FUNC(cgetrs)
#define cheevd_ BLAS_FUNC(cheevd)
#define chemv_ BLAS_FUNC(chemv)
#define cher2_ BLAS_FUNC(cher2)
#define cher2k_ BLAS_FUNC(cher2k)
#define cherk_ BLAS_FUNC(cherk)
#define chetd2_ BLAS_FUNC(chetd2)
#define chetrd_ BLAS_FUNC(chetrd)
#define chseqr_ BLAS_FUNC(chseqr)
#define clabrd_ BLAS_FUNC(clabrd)
#define clacgv_ BLAS_FUNC(clacgv)
#define clacp2_ BLAS_FUNC(clacp2)
#define clacpy_ BLAS_FUNC(clacpy)
#define clacrm_ BLAS_FUNC(clacrm)
#define cladiv_ BLAS_FUNC(cladiv)
#define claed0_ BLAS_FUNC(claed0)
#define claed7_ BLAS_FUNC(claed7)
#define claed8_ BLAS_FUNC(claed8)
#define clahqr_ BLAS_FUNC(clahqr)
#define clahr2_ BLAS_FUNC(clahr2)
#define clals0_ BLAS_FUNC(clals0)
#define clalsa_ BLAS_FUNC(clalsa)
#define clalsd_ BLAS_FUNC(clalsd)
#define clange_ BLAS_FUNC(clange)
#define clanhe_ BLAS_FUNC(clanhe)
#define claqr0_ BLAS_FUNC(claqr0)
#define claqr1_ BLAS_FUNC(claqr1)
#define claqr2_ BLAS_FUNC(claqr2)
#define claqr3_ BLAS_FUNC(claqr3)
#define claqr4_ BLAS_FUNC(claqr4)
#define claqr5_ BLAS_FUNC(claqr5)
#define clarcm_ BLAS_FUNC(clarcm)
#define clarf_ BLAS_FUNC(clarf)
#define clarfb_ BLAS_FUNC(clarfb)
#define clarfg_ BLAS_FUNC(clarfg)
#define clarft_ BLAS_FUNC(clarft)
#define clartg_ BLAS_FUNC(clartg)
#define clascl_ BLAS_FUNC(clascl)
#define claset_ BLAS_FUNC(claset)
#define clasr_ BLAS_FUNC(clasr)
#define classq_ BLAS_FUNC(classq)
#define claswp_ BLAS_FUNC(claswp)
#define clatrd_ BLAS_FUNC(clatrd)
#define clatrs_ BLAS_FUNC(clatrs)
#define clauu2_ BLAS_FUNC(clauu2)
#define clauum_ BLAS_FUNC(clauum)
#define cpotf2_ BLAS_FUNC(cpotf2)
#define cpotrf_ BLAS_FUNC(cpotrf)
#define cpotri_ BLAS_FUNC(cpotri)
#define cpotrs_ BLAS_FUNC(cpotrs)
#define crot_ BLAS_FUNC(crot)
#define cscal_ BLAS_FUNC(cscal)
#define csrot_ BLAS_FUNC(csrot)
#define csscal_ BLAS_FUNC(csscal)
#define cstedc_ BLAS_FUNC(cstedc)
#define csteqr_ BLAS_FUNC(csteqr)
// 定义宏，将函数名映射到 BLAS_FUNC 宏定义的函数名
#define cswap_ BLAS_FUNC(cswap)
#define ctrevc_ BLAS_FUNC(ctrevc)
#define ctrexc_ BLAS_FUNC(ctrexc)
#define ctrmm_ BLAS_FUNC(ctrmm)
#define ctrmv_ BLAS_FUNC(ctrmv)
#define ctrsm_ BLAS_FUNC(ctrsm)
#define ctrsv_ BLAS_FUNC(ctrsv)
#define ctrti2_ BLAS_FUNC(ctrti2)
#define ctrtri_ BLAS_FUNC(ctrtri)
#define cung2r_ BLAS_FUNC(cung2r)
#define cungbr_ BLAS_FUNC(cungbr)
#define cunghr_ BLAS_FUNC(cunghr)
#define cungl2_ BLAS_FUNC(cungl2)
#define cunglq_ BLAS_FUNC(cunglq)
#define cungqr_ BLAS_FUNC(cungqr)
#define cunm2l_ BLAS_FUNC(cunm2l)
#define cunm2r_ BLAS_FUNC(cunm2r)
#define cunmbr_ BLAS_FUNC(cunmbr)
#define cunmhr_ BLAS_FUNC(cunmhr)
#define cunml2_ BLAS_FUNC(cunml2)
#define cunmlq_ BLAS_FUNC(cunmlq)
#define cunmql_ BLAS_FUNC(cunmql)
#define cunmqr_ BLAS_FUNC(cunmqr)
#define cunmtr_ BLAS_FUNC(cunmtr)
#define daxpy_ BLAS_FUNC(daxpy)
#define dbdsdc_ BLAS_FUNC(dbdsdc)
#define dbdsqr_ BLAS_FUNC(dbdsqr)
#define dcabs1_ BLAS_FUNC(dcabs1)
#define dcopy_ BLAS_FUNC(dcopy)
#define ddot_ BLAS_FUNC(ddot)
#define dgebak_ BLAS_FUNC(dgebak)
#define dgebal_ BLAS_FUNC(dgebal)
#define dgebd2_ BLAS_FUNC(dgebd2)
#define dgebrd_ BLAS_FUNC(dgebrd)
#define dgeev_ BLAS_FUNC(dgeev)
#define dgehd2_ BLAS_FUNC(dgehd2)
#define dgehrd_ BLAS_FUNC(dgehrd)
#define dgelq2_ BLAS_FUNC(dgelq2)
#define dgelqf_ BLAS_FUNC(dgelqf)
#define dgelsd_ BLAS_FUNC(dgelsd)
#define dgemm_ BLAS_FUNC(dgemm)
#define dgemv_ BLAS_FUNC(dgemv)
#define dgeqr2_ BLAS_FUNC(dgeqr2)
#define dgeqrf_ BLAS_FUNC(dgeqrf)
#define dger_ BLAS_FUNC(dger)
#define dgesdd_ BLAS_FUNC(dgesdd)
#define dgesv_ BLAS_FUNC(dgesv)
#define dgetf2_ BLAS_FUNC(dgetf2)
#define dgetrf_ BLAS_FUNC(dgetrf)
#define dgetrs_ BLAS_FUNC(dgetrs)
#define dhseqr_ BLAS_FUNC(dhseqr)
#define disnan_ BLAS_FUNC(disnan)
#define dlabad_ BLAS_FUNC(dlabad)
#define dlabrd_ BLAS_FUNC(dlabrd)
#define dlacpy_ BLAS_FUNC(dlacpy)
#define dladiv_ BLAS_FUNC(dladiv)
#define dlae2_ BLAS_FUNC(dlae2)
#define dlaed0_ BLAS_FUNC(dlaed0)
#define dlaed1_ BLAS_FUNC(dlaed1)
#define dlaed2_ BLAS_FUNC(dlaed2)
#define dlaed3_ BLAS_FUNC(dlaed3)
#define dlaed4_ BLAS_FUNC(dlaed4)
#define dlaed5_ BLAS_FUNC(dlaed5)
#define dlaed6_ BLAS_FUNC(dlaed6)
#define dlaed7_ BLAS_FUNC(dlaed7)
#define dlaed8_ BLAS_FUNC(dlaed8)
#define dlaed9_ BLAS_FUNC(dlaed9)
#define dlaeda_ BLAS_FUNC(dlaeda)
#define dlaev2_ BLAS_FUNC(dlaev2)
#define dlaexc_ BLAS_FUNC(dlaexc)
#define dlahqr_ BLAS_FUNC(dlahqr)
#define dlahr2_ BLAS_FUNC(dlahr2)
#define dlaisnan_ BLAS_FUNC(dlaisnan)
#define dlaln2_ BLAS_FUNC(dlaln2)
#define dlals0_ BLAS_FUNC(dlals0)
#define dlalsa_ BLAS_FUNC(dlalsa)
#define dlalsd_ BLAS_FUNC(dlalsd)
#define dlamc1_ BLAS_FUNC(dlamc1)
#define dlamc2_ BLAS_FUNC(dlamc2)
#define dlamc3_ BLAS_FUNC(dlamc3)
#define dlamc4_ BLAS_FUNC(dlamc4)
#define dlamc5_ BLAS_FUNC(dlamc5)
#define dlamch_ BLAS_FUNC(dlamch)
#define dlamrg_ BLAS_FUNC(dlamrg)
#define dlange_ BLAS_FUNC(dlange)
#define dlanst_ BLAS_FUNC(dlanst)
#define dlansy_ BLAS_FUNC(dlansy)
#define dlanv2_ BLAS_FUNC(dlanv2)
#define dlapy2_ BLAS_FUNC(dlapy2)
// 定义各个 BLAS 函数的宏，它们是对应于 LAPACK 和 BLAS 数学库中特定函数的符号名称映射
#define dlapy3_ BLAS_FUNC(dlapy3)
#define dlaqr0_ BLAS_FUNC(dlaqr0)
#define dlaqr1_ BLAS_FUNC(dlaqr1)
#define dlaqr2_ BLAS_FUNC(dlaqr2)
#define dlaqr3_ BLAS_FUNC(dlaqr3)
#define dlaqr4_ BLAS_FUNC(dlaqr4)
#define dlaqr5_ BLAS_FUNC(dlaqr5)
#define dlarf_ BLAS_FUNC(dlarf)
#define dlarfb_ BLAS_FUNC(dlarfb)
#define dlarfg_ BLAS_FUNC(dlarfg)
#define dlarft_ BLAS_FUNC(dlarft)
#define dlarfx_ BLAS_FUNC(dlarfx)
#define dlartg_ BLAS_FUNC(dlartg)
#define dlas2_ BLAS_FUNC(dlas2)
#define dlascl_ BLAS_FUNC(dlascl)
#define dlasd0_ BLAS_FUNC(dlasd0)
#define dlasd1_ BLAS_FUNC(dlasd1)
#define dlasd2_ BLAS_FUNC(dlasd2)
#define dlasd3_ BLAS_FUNC(dlasd3)
#define dlasd4_ BLAS_FUNC(dlasd4)
#define dlasd5_ BLAS_FUNC(dlasd5)
#define dlasd6_ BLAS_FUNC(dlasd6)
#define dlasd7_ BLAS_FUNC(dlasd7)
#define dlasd8_ BLAS_FUNC(dlasd8)
#define dlasda_ BLAS_FUNC(dlasda)
#define dlasdq_ BLAS_FUNC(dlasdq)
#define dlasdt_ BLAS_FUNC(dlasdt)
#define dlaset_ BLAS_FUNC(dlaset)
#define dlasq1_ BLAS_FUNC(dlasq1)
#define dlasq2_ BLAS_FUNC(dlasq2)
#define dlasq3_ BLAS_FUNC(dlasq3)
#define dlasq4_ BLAS_FUNC(dlasq4)
#define dlasq5_ BLAS_FUNC(dlasq5)
#define dlasq6_ BLAS_FUNC(dlasq6)
#define dlasr_ BLAS_FUNC(dlasr)
#define dlasrt_ BLAS_FUNC(dlasrt)
#define dlassq_ BLAS_FUNC(dlassq)
#define dlasv2_ BLAS_FUNC(dlasv2)
#define dlaswp_ BLAS_FUNC(dlaswp)
#define dlasy2_ BLAS_FUNC(dlasy2)
#define dlatrd_ BLAS_FUNC(dlatrd)
#define dlauu2_ BLAS_FUNC(dlauu2)
#define dlauum_ BLAS_FUNC(dlauum)
#define dnrm2_ BLAS_FUNC(dnrm2)
#define dorg2r_ BLAS_FUNC(dorg2r)
#define dorgbr_ BLAS_FUNC(dorgbr)
#define dorghr_ BLAS_FUNC(dorghr)
#define dorgl2_ BLAS_FUNC(dorgl2)
#define dorglq_ BLAS_FUNC(dorglq)
#define dorgqr_ BLAS_FUNC(dorgqr)
#define dorm2l_ BLAS_FUNC(dorm2l)
#define dorm2r_ BLAS_FUNC(dorm2r)
#define dormbr_ BLAS_FUNC(dormbr)
#define dormhr_ BLAS_FUNC(dormhr)
#define dorml2_ BLAS_FUNC(dorml2)
#define dormlq_ BLAS_FUNC(dormlq)
#define dormql_ BLAS_FUNC(dormql)
#define dormqr_ BLAS_FUNC(dormqr)
#define dormtr_ BLAS_FUNC(dormtr)
#define dpotf2_ BLAS_FUNC(dpotf2)
#define dpotrf_ BLAS_FUNC(dpotrf)
#define dpotri_ BLAS_FUNC(dpotri)
#define dpotrs_ BLAS_FUNC(dpotrs)
#define drot_ BLAS_FUNC(drot)
#define dscal_ BLAS_FUNC(dscal)
#define dstedc_ BLAS_FUNC(dstedc)
#define dsteqr_ BLAS_FUNC(dsteqr)
#define dsterf_ BLAS_FUNC(dsterf)
#define dswap_ BLAS_FUNC(dswap)
#define dsyevd_ BLAS_FUNC(dsyevd)
#define dsymv_ BLAS_FUNC(dsymv)
#define dsyr2_ BLAS_FUNC(dsyr2)
#define dsyr2k_ BLAS_FUNC(dsyr2k)
#define dsyrk_ BLAS_FUNC(dsyrk)
#define dsytd2_ BLAS_FUNC(dsytd2)
#define dsytrd_ BLAS_FUNC(dsytrd)
#define dtrevc_ BLAS_FUNC(dtrevc)
#define dtrexc_ BLAS_FUNC(dtrexc)
#define dtrmm_ BLAS_FUNC(dtrmm)
#define dtrmv_ BLAS_FUNC(dtrmv)
#define dtrsm_ BLAS_FUNC(dtrsm)
#define dtrti2_ BLAS_FUNC(dtrti2)
#define dtrtri_ BLAS_FUNC(dtrtri)
#define dzasum_ BLAS_FUNC(dzasum)
#define dznrm2_ BLAS_FUNC(dznrm2)
#define icamax_ BLAS_FUNC(icamax)
#define idamax_ BLAS_FUNC(idamax)
#define ieeeck_ BLAS_FUNC(ieeeck)
#define ilaclc_ BLAS_FUNC(ilaclc)
#define ilaclr_ BLAS_FUNC(ilaclr)
// 定义宏 ilaclr_，展开为 BLAS_FUNC(ilaclr) 的内容
#define iladlc_ BLAS_FUNC(iladlc)
// 定义宏 iladlc_，展开为 BLAS_FUNC(iladlc) 的内容
#define iladlr_ BLAS_FUNC(iladlr)
// 定义宏 iladlr_，展开为 BLAS_FUNC(iladlr) 的内容
#define ilaenv_ BLAS_FUNC(ilaenv)
// 定义宏 ilaenv_，展开为 BLAS_FUNC(ilaenv) 的内容
#define ilaslc_ BLAS_FUNC(ilaslc)
// 定义宏 ilaslc_，展开为 BLAS_FUNC(ilaslc) 的内容
#define ilaslr_ BLAS_FUNC(ilaslr)
// 定义宏 ilaslr_，展开为 BLAS_FUNC(ilaslr) 的内容
#define ilazlc_ BLAS_FUNC(ilazlc)
// 定义宏 ilazlc_，展开为 BLAS_FUNC(ilazlc) 的内容
#define ilazlr_ BLAS_FUNC(ilazlr)
// 定义宏 ilazlr_，展开为 BLAS_FUNC(ilazlr) 的内容
#define iparmq_ BLAS_FUNC(iparmq)
// 定义宏 iparmq_，展开为 BLAS_FUNC(iparmq) 的内容
#define isamax_ BLAS_FUNC(isamax)
// 定义宏 isamax_，展开为 BLAS_FUNC(isamax) 的内容
#define izamax_ BLAS_FUNC(izamax)
// 定义宏 izamax_，展开为 BLAS_FUNC(izamax) 的内容
#define lsame_ BLAS_FUNC(lsame)
// 定义宏 lsame_，展开为 BLAS_FUNC(lsame) 的内容
#define saxpy_ BLAS_FUNC(saxpy)
// 定义宏 saxpy_，展开为 BLAS_FUNC(saxpy) 的内容
#define sbdsdc_ BLAS_FUNC(sbdsdc)
// 定义宏 sbdsdc_，展开为 BLAS_FUNC(sbdsdc) 的内容
#define sbdsqr_ BLAS_FUNC(sbdsqr)
// 定义宏 sbdsqr_，展开为 BLAS_FUNC(sbdsqr) 的内容
#define scabs1_ BLAS_FUNC(scabs1)
// 定义宏 scabs1_，展开为 BLAS_FUNC(scabs1) 的内容
#define scasum_ BLAS_FUNC(scasum)
// 定义宏 scasum_，展开为 BLAS_FUNC(scasum) 的内容
#define scnrm2_ BLAS_FUNC(scnrm2)
// 定义宏 scnrm2_，展开为 BLAS_FUNC(scnrm2) 的内容
#define scopy_ BLAS_FUNC(scopy)
// 定义宏 scopy_，展开为 BLAS_FUNC(scopy) 的内容
#define sdot_ BLAS_FUNC(sdot)
// 定义宏 sdot_，展开为 BLAS_FUNC(sdot) 的内容
#define sgebak_ BLAS_FUNC(sgebak)
// 定义宏 sgebak_，展开为 BLAS_FUNC(sgebak) 的内容
#define sgebal_ BLAS_FUNC(sgebal)
// 定义宏 sgebal_，展开为 BLAS_FUNC(sgebal) 的内容
#define sgebd2_ BLAS_FUNC(sgebd2)
// 定义宏 sgebd2_，展开为 BLAS_FUNC(sgebd2) 的内容
#define sgebrd_ BLAS_FUNC(sgebrd)
// 定义宏 sgebrd_，展开为 BLAS_FUNC(sgebrd) 的内容
#define sgeev_ BLAS_FUNC(sgeev)
// 定义宏 sgeev_，展开为 BLAS_FUNC(sgeev) 的内容
#define sgehd2_ BLAS_FUNC(sgehd2)
// 定义宏 sgehd2_，展开为 BLAS_FUNC(sgehd2) 的内容
#define sgehrd_ BLAS_FUNC(sgehrd)
// 定义宏 sgehrd_，展开为 BLAS_FUNC(sgehrd) 的内容
#define sgelq2_ BLAS_FUNC(sgelq2)
// 定义宏 sgelq2_，展开为 BLAS_FUNC(sgelq2) 的内容
#define sgelqf_ BLAS_FUNC(sgelqf)
// 定义宏 sgelqf_，展开为 BLAS_FUNC(sgelqf) 的内容
#define sgelsd_ BLAS_FUNC(sgelsd)
// 定义宏 sgelsd_，展开为 BLAS_FUNC(sgelsd) 的内容
#define sgemm_ BLAS_FUNC(sgemm)
// 定义宏 sgemm_，展开为 BLAS_FUNC(sgemm) 的内容
#define sgemv_ BLAS_FUNC(sgemv)
// 定义宏 sgemv_，展开为 BLAS_FUNC(sgemv) 的内容
#define sgeqr2_ BLAS_FUNC(sgeqr2)
// 定义宏 sgeqr2_，展开为 BLAS_FUNC(sgeqr2) 的内容
#define sgeqrf_ BLAS_FUNC(sgeqrf)
// 定义宏 sgeqrf_，展开为 BLAS_FUNC(sgeqrf) 的内容
#define sger_ BLAS_FUNC(sger)
// 定义宏 sger_，展开为 BLAS_FUNC(sger) 的内容
#define sgesdd_ BLAS_FUNC(sgesdd)
// 定义宏 sgesdd_，展开为 BLAS_FUNC(sgesdd) 的内容
#define sgesv_ BLAS_FUNC(sgesv)
// 定义宏 sgesv_，展开为 BLAS_FUNC(sgesv) 的内容
#define sgetf2_ BLAS_FUNC(sgetf2)
// 定义宏 sgetf2_，展开为 BLAS_FUNC(sgetf2) 的内容
#define sgetrf_ BLAS_FUNC(sgetrf)
// 定义宏 sgetrf_，展开为 BLAS_FUNC(sgetrf) 的内容
#define sgetrs_ BLAS_FUNC(sgetrs)
// 定义宏 sgetrs_，展开为 BLAS_FUNC(sgetrs) 的内容
#define shseqr_ BLAS_FUNC(shseqr)
// 定义宏 shseqr_，展开为 BLAS_FUNC(shseqr) 的内容
#define sisnan_ BLAS_FUNC(sisnan)
// 定义宏 sisnan_，展开为 BLAS_FUNC(sisnan) 的内容
#define slabad_ BLAS_FUNC(slabad)
// 定义宏 slabad_，展开为 BLAS_FUNC(slabad) 的内容
#define slabrd_ BLAS_FUNC(slabrd)
// 定义宏 slabrd_，展开为 BLAS_FUNC(slabrd) 的内容
#define slacpy_ BLAS_FUNC(slacpy)
// 定义宏 slacpy_，展开为 BLAS_FUNC(slacpy) 的内容
#define sladiv_ BLAS_FUNC(sladiv)
// 定义宏 sladiv_，展开为 BLAS_FUNC(sladiv) 的内容
#define slae2_ BLAS_FUNC(slae2)
// 定义宏 slae2_，展开为 BLAS_FUNC(slae2) 的内容
#define slaed0_ BLAS_FUNC(slaed0)
// 定义宏 slaed0_，展开为 BLAS_FUNC(slaed0) 的内容
#define slaed1_ BLAS_FUNC(slaed1)
// 定义宏 slaed1_，展开为 BLAS_FUNC(slaed1) 的内容
#define slaed2_ BLAS_FUNC(slaed2)
// 定义宏 slaed2_，展开为 BLAS_FUNC(slaed2) 的内容
#define slaed3_ BLAS_FUNC(slaed3)
// 定义宏 slaed3_，展开为 BLAS_FUNC(slaed3) 的内容
#define slaed4_ BLAS_FUNC(slaed4)
// 定义宏 slaed4_，展开为 BLAS_FUNC(slaed4) 的内容
#define slaed5_ BLAS_FUNC(slaed5)
// 定义宏 slaed5_，展开为 BLAS_FUNC(slaed5) 的内容
#define slaed6_ BLAS_FUNC(slaed6)
// 定义宏 slaed6_，展开为 BLAS_FUNC(slaed6) 的内容
#define slaed7_ BLAS_FUNC(slaed7)
// 定义宏 slaed7_，展开为 BLAS_FUNC(slaed7) 的内容
#define slaed8_ BLAS_FUNC(slaed8)
// 定义宏 slaed8_，展开为 BLAS_FUNC(slaed8) 的内容
#define slaed
// 定义各个 BLAS 函数的别名，使用 BLAS_FUNC 宏进行定义
#define slarfg_ BLAS_FUNC(slarfg)
#define slarft_ BLAS_FUNC(slarft)
#define slarfx_ BLAS_FUNC(slarfx)
#define slartg_ BLAS_FUNC(slartg)
#define slas2_ BLAS_FUNC(slas2)
#define slascl_ BLAS_FUNC(slascl)
#define slasd0_ BLAS_FUNC(slasd0)
#define slasd1_ BLAS_FUNC(slasd1)
#define slasd2_ BLAS_FUNC(slasd2)
#define slasd3_ BLAS_FUNC(slasd3)
#define slasd4_ BLAS_FUNC(slasd4)
#define slasd5_ BLAS_FUNC(slasd5)
#define slasd6_ BLAS_FUNC(slasd6)
#define slasd7_ BLAS_FUNC(slasd7)
#define slasd8_ BLAS_FUNC(slasd8)
#define slasda_ BLAS_FUNC(slasda)
#define slasdq_ BLAS_FUNC(slasdq)
#define slasdt_ BLAS_FUNC(slasdt)
#define slaset_ BLAS_FUNC(slaset)
#define slasq1_ BLAS_FUNC(slasq1)
#define slasq2_ BLAS_FUNC(slasq2)
#define slasq3_ BLAS_FUNC(slasq3)
#define slasq4_ BLAS_FUNC(slasq4)
#define slasq5_ BLAS_FUNC(slasq5)
#define slasq6_ BLAS_FUNC(slasq6)
#define slasr_ BLAS_FUNC(slasr)
#define slasrt_ BLAS_FUNC(slasrt)
#define slassq_ BLAS_FUNC(slassq)
#define slasv2_ BLAS_FUNC(slasv2)
#define slaswp_ BLAS_FUNC(slaswp)
#define slasy2_ BLAS_FUNC(slasy2)
#define slatrd_ BLAS_FUNC(slatrd)
#define slauu2_ BLAS_FUNC(slauu2)
#define slauum_ BLAS_FUNC(slauum)
#define snrm2_ BLAS_FUNC(snrm2)
#define sorg2r_ BLAS_FUNC(sorg2r)
#define sorgbr_ BLAS_FUNC(sorgbr)
#define sorghr_ BLAS_FUNC(sorghr)
#define sorgl2_ BLAS_FUNC(sorgl2)
#define sorglq_ BLAS_FUNC(sorglq)
#define sorgqr_ BLAS_FUNC(sorgqr)
#define sorm2l_ BLAS_FUNC(sorm2l)
#define sorm2r_ BLAS_FUNC(sorm2r)
#define sormbr_ BLAS_FUNC(sormbr)
#define sormhr_ BLAS_FUNC(sormhr)
#define sorml2_ BLAS_FUNC(sorml2)
#define sormlq_ BLAS_FUNC(sormlq)
#define sormql_ BLAS_FUNC(sormql)
#define sormqr_ BLAS_FUNC(sormqr)
#define sormtr_ BLAS_FUNC(sormtr)
#define spotf2_ BLAS_FUNC(spotf2)
#define spotrf_ BLAS_FUNC(spotrf)
#define spotri_ BLAS_FUNC(spotri)
#define spotrs_ BLAS_FUNC(spotrs)
#define srot_ BLAS_FUNC(srot)
#define sscal_ BLAS_FUNC(sscal)
#define sstedc_ BLAS_FUNC(sstedc)
#define ssteqr_ BLAS_FUNC(ssteqr)
#define ssterf_ BLAS_FUNC(ssterf)
#define sswap_ BLAS_FUNC(sswap)
#define ssyevd_ BLAS_FUNC(ssyevd)
#define ssymv_ BLAS_FUNC(ssymv)
#define ssyr2_ BLAS_FUNC(ssyr2)
#define ssyr2k_ BLAS_FUNC(ssyr2k)
#define ssyrk_ BLAS_FUNC(ssyrk)
#define ssytd2_ BLAS_FUNC(ssytd2)
#define ssytrd_ BLAS_FUNC(ssytrd)
#define strevc_ BLAS_FUNC(strevc)
#define strexc_ BLAS_FUNC(strexc)
#define strmm_ BLAS_FUNC(strmm)
#define strmv_ BLAS_FUNC(strmv)
#define strsm_ BLAS_FUNC(strsm)
#define strti2_ BLAS_FUNC(strti2)
#define strtri_ BLAS_FUNC(strtri)
#define xerbla_ BLAS_FUNC(xerbla)
#define zaxpy_ BLAS_FUNC(zaxpy)
#define zcopy_ BLAS_FUNC(zcopy)
#define zdotc_ BLAS_FUNC(zdotc)
#define zdotu_ BLAS_FUNC(zdotu)
#define zdrot_ BLAS_FUNC(zdrot)
#define zdscal_ BLAS_FUNC(zdscal)
#define zgebak_ BLAS_FUNC(zgebak)
#define zgebal_ BLAS_FUNC(zgebal)
#define zgebd2_ BLAS_FUNC(zgebd2)
#define zgebrd_ BLAS_FUNC(zgebrd)
#define zgeev_ BLAS_FUNC(zgeev)
#define zgehd2_ BLAS_FUNC(zgehd2)
#define zgehrd_ BLAS_FUNC(zgehrd)
#define zgelq2_ BLAS_FUNC(zgelq2)
// 定义宏，用于重命名 BLAS 函数名，以下是一系列宏定义
#define zgelqf_ BLAS_FUNC(zgelqf)
#define zgelsd_ BLAS_FUNC(zgelsd)
#define zgemm_ BLAS_FUNC(zgemm)
#define zgemv_ BLAS_FUNC(zgemv)
#define zgeqr2_ BLAS_FUNC(zgeqr2)
#define zgeqrf_ BLAS_FUNC(zgeqrf)
#define zgerc_ BLAS_FUNC(zgerc)
#define zgeru_ BLAS_FUNC(zgeru)
#define zgesdd_ BLAS_FUNC(zgesdd)
#define zgesv_ BLAS_FUNC(zgesv)
#define zgetf2_ BLAS_FUNC(zgetf2)
#define zgetrf_ BLAS_FUNC(zgetrf)
#define zgetrs_ BLAS_FUNC(zgetrs)
#define zheevd_ BLAS_FUNC(zheevd)
#define zhemv_ BLAS_FUNC(zhemv)
#define zher2_ BLAS_FUNC(zher2)
#define zher2k_ BLAS_FUNC(zher2k)
#define zherk_ BLAS_FUNC(zherk)
#define zhetd2_ BLAS_FUNC(zhetd2)
#define zhetrd_ BLAS_FUNC(zhetrd)
#define zhseqr_ BLAS_FUNC(zhseqr)
#define zlabrd_ BLAS_FUNC(zlabrd)
#define zlacgv_ BLAS_FUNC(zlacgv)
#define zlacp2_ BLAS_FUNC(zlacp2)
#define zlacpy_ BLAS_FUNC(zlacpy)
#define zlacrm_ BLAS_FUNC(zlacrm)
#define zladiv_ BLAS_FUNC(zladiv)
#define zlaed0_ BLAS_FUNC(zlaed0)
#define zlaed7_ BLAS_FUNC(zlaed7)
#define zlaed8_ BLAS_FUNC(zlaed8)
#define zlahqr_ BLAS_FUNC(zlahqr)
#define zlahr2_ BLAS_FUNC(zlahr2)
#define zlals0_ BLAS_FUNC(zlals0)
#define zlalsa_ BLAS_FUNC(zlalsa)
#define zlalsd_ BLAS_FUNC(zlalsd)
#define zlange_ BLAS_FUNC(zlange)
#define zlanhe_ BLAS_FUNC(zlanhe)
#define zlaqr0_ BLAS_FUNC(zlaqr0)
#define zlaqr1_ BLAS_FUNC(zlaqr1)
#define zlaqr2_ BLAS_FUNC(zlaqr2)
#define zlaqr3_ BLAS_FUNC(zlaqr3)
#define zlaqr4_ BLAS_FUNC(zlaqr4)
#define zlaqr5_ BLAS_FUNC(zlaqr5)
#define zlarcm_ BLAS_FUNC(zlarcm)
#define zlarf_ BLAS_FUNC(zlarf)
#define zlarfb_ BLAS_FUNC(zlarfb)
#define zlarfg_ BLAS_FUNC(zlarfg)
#define zlarft_ BLAS_FUNC(zlarft)
#define zlartg_ BLAS_FUNC(zlartg)
#define zlascl_ BLAS_FUNC(zlascl)
#define zlaset_ BLAS_FUNC(zlaset)
#define zlasr_ BLAS_FUNC(zlasr)
#define zlassq_ BLAS_FUNC(zlassq)
#define zlaswp_ BLAS_FUNC(zlaswp)
#define zlatrd_ BLAS_FUNC(zlatrd)
#define zlatrs_ BLAS_FUNC(zlatrs)
#define zlauu2_ BLAS_FUNC(zlauu2)
#define zlauum_ BLAS_FUNC(zlauum)
#define zpotf2_ BLAS_FUNC(zpotf2)
#define zpotrf_ BLAS_FUNC(zpotrf)
#define zpotri_ BLAS_FUNC(zpotri)
#define zpotrs_ BLAS_FUNC(zpotrs)
#define zrot_ BLAS_FUNC(zrot)
#define zscal_ BLAS_FUNC(zscal)
#define zstedc_ BLAS_FUNC(zstedc)
#define zsteqr_ BLAS_FUNC(zsteqr)
#define zswap_ BLAS_FUNC(zswap)
#define ztrevc_ BLAS_FUNC(ztrevc)
#define ztrexc_ BLAS_FUNC(ztrexc)
#define ztrmm_ BLAS_FUNC(ztrmm)
#define ztrmv_ BLAS_FUNC(ztrmv)
#define ztrsm_ BLAS_FUNC(ztrsm)
#define ztrsv_ BLAS_FUNC(ztrsv)
#define ztrti2_ BLAS_FUNC(ztrti2)
#define ztrtri_ BLAS_FUNC(ztrtri)
#define zung2r_ BLAS_FUNC(zung2r)
#define zungbr_ BLAS_FUNC(zungbr)
#define zunghr_ BLAS_FUNC(zunghr)
#define zungl2_ BLAS_FUNC(zungl2)
#define zunglq_ BLAS_FUNC(zunglq)
#define zungqr_ BLAS_FUNC(zungqr)
#define zunm2l_ BLAS_FUNC(zunm2l)
#define zunm2r_ BLAS_FUNC(zunm2r)
#define zunmbr_ BLAS_FUNC(zunmbr)
#define zunmhr_ BLAS_FUNC(zunmhr)
#define zunml2_ BLAS_FUNC(zunml2)
#define zunmlq_ BLAS_FUNC(zunmlq)
#define zunmql_ BLAS_FUNC(zunmql)
#define zunmqr_ BLAS_FUNC(zunmqr)
# 定义 zunmtr_ 符号，表示 BLAS_FUNC 函数的 zunmtr 实现
#define zunmtr_ BLAS_FUNC(zunmtr)
# 下面是一系列由 f2c.c 导出的符号名称重定义，这些符号名称都指向 numpy_lapack_lite 模块中对应的函数或变量。
# 例如，将 abort_ 重定义为 numpy_lapack_lite_abort_
#define abort_ numpy_lapack_lite_abort_
#define c_abs numpy_lapack_lite_c_abs
#define c_cos numpy_lapack_lite_c_cos
#define c_div numpy_lapack_lite_c_div
#define c_exp numpy_lapack_lite_c_exp
#define c_log numpy_lapack_lite_c_log
#define c_sin numpy_lapack_lite_c_sin
#define c_sqrt numpy_lapack_lite_c_sqrt
#define d_abs numpy_lapack_lite_d_abs
#define d_acos numpy_lapack_lite_d_acos
#define d_asin numpy_lapack_lite_d_asin
#define d_atan numpy_lapack_lite_d_atan
#define d_atn2 numpy_lapack_lite_d_atn2
#define d_cnjg numpy_lapack_lite_d_cnjg
#define d_cos numpy_lapack_lite_d_cos
#define d_cosh numpy_lapack_lite_d_cosh
#define d_dim numpy_lapack_lite_d_dim
#define d_exp numpy_lapack_lite_d_exp
#define d_imag numpy_lapack_lite_d_imag
#define d_int numpy_lapack_lite_d_int
#define d_lg10 numpy_lapack_lite_d_lg10
#define d_log numpy_lapack_lite_d_log
#define d_mod numpy_lapack_lite_d_mod
#define d_nint numpy_lapack_lite_d_nint
#define d_prod numpy_lapack_lite_d_prod
#define d_sign numpy_lapack_lite_d_sign
#define d_sin numpy_lapack_lite_d_sin
#define d_sinh numpy_lapack_lite_d_sinh
#define d_sqrt numpy_lapack_lite_d_sqrt
#define d_tan numpy_lapack_lite_d_tan
#define d_tanh numpy_lapack_lite_d_tanh
#define derf_ numpy_lapack_lite_derf_
#define derfc_ numpy_lapack_lite_derfc_
#define do_fio numpy_lapack_lite_do_fio
#define do_lio numpy_lapack_lite_do_lio
#define do_uio numpy_lapack_lite_do_uio
#define e_rdfe numpy_lapack_lite_e_rdfe
#define e_rdue numpy_lapack_lite_e_rdue
#define e_rsfe numpy_lapack_lite_e_rsfe
#define e_rsfi numpy_lapack_lite_e_rsfi
#define e_rsle numpy_lapack_lite_e_rsle
#define e_rsli numpy_lapack_lite_e_rsli
#define e_rsue numpy_lapack_lite_e_rsue
#define e_wdfe numpy_lapack_lite_e_wdfe
#define e_wdue numpy_lapack_lite_e_wdue
#define e_wsfe numpy_lapack_lite_e_wsfe
#define e_wsfi numpy_lapack_lite_e_wsfi
#define e_wsle numpy_lapack_lite_e_wsle
#define e_wsli numpy_lapack_lite_e_wsli
#define e_wsue numpy_lapack_lite_e_wsue
#define ef1asc_ numpy_lapack_lite_ef1asc_
#define ef1cmc_ numpy_lapack_lite_ef1cmc_
#define erf_ numpy_lapack_lite_erf_
#define erfc_ numpy_lapack_lite_erfc_
#define f__cabs numpy_lapack_lite_f__cabs
#define f__cabsf numpy_lapack_lite_f__cabsf
#define f_back numpy_lapack_lite_f_back
#define f_clos numpy_lapack_lite_f_clos
#define f_end numpy_lapack_lite_f_end
#define f_exit numpy_lapack_lite_f_exit
#define f_inqu numpy_lapack_lite_f_inqu
#define f_open numpy_lapack_lite_f_open
#define f_rew numpy_lapack_lite_f_rew
#define flush_ numpy_lapack_lite_flush_
#define getarg_ numpy_lapack_lite_getarg_
#define getenv_ numpy_lapack_lite_getenv_
#define h_abs numpy_lapack_lite_h_abs
#define h_dim numpy_lapack_lite_h_dim
#define h_dnnt numpy_lapack_lite_h_dnnt
#define h_indx numpy_lapack_lite_h_indx
#define h_len numpy_lapack_lite_h_len
#define h_mod numpy_lapack_lite_h_mod
#define h_nint numpy_lapack_lite_h_nint
#define h_sign numpy_lapack_lite_h_sign
# 定义宏指令，用于将给定名称重新映射为对应的numpy_lapack_lite模块中的函数或变量
#define hl_ge numpy_lapack_lite_hl_ge
#define hl_gt numpy_lapack_lite_hl_gt
#define hl_le numpy_lapack_lite_hl_le
#define hl_lt numpy_lapack_lite_hl_lt
#define i_abs numpy_lapack_lite_i_abs
#define i_dim numpy_lapack_lite_i_dim
#define i_dnnt numpy_lapack_lite_i_dnnt
#define i_indx numpy_lapack_lite_i_indx
#define i_len numpy_lapack_lite_i_len
#define i_mod numpy_lapack_lite_i_mod
#define i_nint numpy_lapack_lite_i_nint
#define i_sign numpy_lapack_lite_i_sign
#define iargc_ numpy_lapack_lite_iargc_
#define l_ge numpy_lapack_lite_l_ge
#define l_gt numpy_lapack_lite_l_gt
#define l_le numpy_lapack_lite_l_le
#define l_lt numpy_lapack_lite_l_lt
#define pow_ci numpy_lapack_lite_pow_ci
#define pow_dd numpy_lapack_lite_pow_dd
#define pow_di numpy_lapack_lite_pow_di
#define pow_hh numpy_lapack_lite_pow_hh
#define pow_ii numpy_lapack_lite_pow_ii
#define pow_ri numpy_lapack_lite_pow_ri
#define pow_zi numpy_lapack_lite_pow_zi
#define pow_zz numpy_lapack_lite_pow_zz
#define r_abs numpy_lapack_lite_r_abs
#define r_acos numpy_lapack_lite_r_acos
#define r_asin numpy_lapack_lite_r_asin
#define r_atan numpy_lapack_lite_r_atan
#define r_atn2 numpy_lapack_lite_r_atn2
#define r_cnjg numpy_lapack_lite_r_cnjg
#define r_cos numpy_lapack_lite_r_cos
#define r_cosh numpy_lapack_lite_r_cosh
#define r_dim numpy_lapack_lite_r_dim
#define r_exp numpy_lapack_lite_r_exp
#define r_imag numpy_lapack_lite_r_imag
#define r_int numpy_lapack_lite_r_int
#define r_lg10 numpy_lapack_lite_r_lg10
#define r_log numpy_lapack_lite_r_log
#define r_mod numpy_lapack_lite_r_mod
#define r_nint numpy_lapack_lite_r_nint
#define r_sign numpy_lapack_lite_r_sign
#define r_sin numpy_lapack_lite_r_sin
#define r_sinh numpy_lapack_lite_r_sinh
#define r_sqrt numpy_lapack_lite_r_sqrt
#define r_tan numpy_lapack_lite_r_tan
#define r_tanh numpy_lapack_lite_r_tanh
#define s_cat numpy_lapack_lite_s_cat
#define s_cmp numpy_lapack_lite_s_cmp
#define s_copy numpy_lapack_lite_s_copy
#define s_paus numpy_lapack_lite_s_paus
#define s_rdfe numpy_lapack_lite_s_rdfe
#define s_rdue numpy_lapack_lite_s_rdue
#define s_rnge numpy_lapack_lite_s_rnge
#define s_rsfe numpy_lapack_lite_s_rsfe
#define s_rsfi numpy_lapack_lite_s_rsfi
#define s_rsle numpy_lapack_lite_s_rsle
#define s_rsli numpy_lapack_lite_s_rsli
#define s_rsne numpy_lapack_lite_s_rsne
#define s_rsni numpy_lapack_lite_s_rsni
#define s_rsue numpy_lapack_lite_s_rsue
#define s_stop numpy_lapack_lite_s_stop
#define s_wdfe numpy_lapack_lite_s_wdfe
#define s_wdue numpy_lapack_lite_s_wdue
#define s_wsfe numpy_lapack_lite_s_wsfe
#define s_wsfi numpy_lapack_lite_s_wsfi
#define s_wsle numpy_lapack_lite_s_wsle
#define s_wsli numpy_lapack_lite_s_wsli
#define s_wsne numpy_lapack_lite_s_wsne
#define s_wsni numpy_lapack_lite_s_wsni
#define s_wsue numpy_lapack_lite_s_wsue
#define sig_die numpy_lapack_lite_sig_die
#define signal_ numpy_lapack_lite_signal_
#define system_ numpy_lapack_lite_system_
#define z_abs numpy_lapack_lite_z_abs
#define z_cos numpy_lapack_lite_z_cos
# 定义五个宏，分别将其映射到 numpy_lapack_lite 库中的对应函数
#define z_div numpy_lapack_lite_z_div
#define z_exp numpy_lapack_lite_z_exp
#define z_log numpy_lapack_lite_z_log
#define z_sin numpy_lapack_lite_z_sin
#define z_sqrt numpy_lapack_lite_z_sqrt
.\numpy\numpy\linalg\lapack_lite\make_lite.py
#!/usr/bin/env python2.7
# WARNING! This a Python 2 script. Read README.rst for rationale.
"""
Usage: make_lite.py <wrapped_routines_file> <lapack_dir>
Typical invocation:
make_lite.py wrapped_routines /tmp/lapack-3.x.x
Requires the following to be on the path:
* f2c
* patch
"""
import sys # 导入 sys 模块，用于访问命令行参数和系统相关功能
import os # 导入 os 模块，提供了对操作系统进行调用的接口
import re # 导入 re 模块，用于处理正则表达式的匹配和操作
import subprocess # 导入 subprocess 模块，用于执行外部命令
import shutil # 导入 shutil 模块，提供了一些高级的文件操作功能
import fortran # 导入自定义的 fortran 模块，用于处理 Fortran 相关操作
import clapack_scrub # 导入自定义的 clapack_scrub 模块，用于处理 Clapack 相关操作
try:
from distutils.spawn import find_executable as which # 尝试从 distutils.spawn 模块导入 find_executable 函数（Python 2）
except ImportError:
from shutil import which # 如果导入失败，则从 shutil 模块导入 which 函数（Python 3）
# Arguments to pass to f2c. You'll always want -A for ANSI C prototypes
# Others of interest: -a to not make variables static by default
# -C to check array subscripts
F2C_ARGS = ['-A', '-Nx800'] # 设置传递给 f2c 的参数列表
# The header to add to the top of the f2c_*.c file. Note that dlamch_() calls
# will be replaced by the macros below by clapack_scrub.scrub_source()
HEADER_BLURB = '''\
/*
* NOTE: This is generated code. Look in numpy/linalg/lapack_lite for
* information on remaking this file.
*/
'''
HEADER = HEADER_BLURB + '''\
#include "f2c.h"
#ifdef HAVE_CONFIG
#include "config.h"
#else
extern doublereal dlamch_(char *);
#define EPSILON dlamch_("Epsilon")
#define SAFEMINIMUM dlamch_("Safe minimum")
#define PRECISION dlamch_("Precision")
#define BASE dlamch_("Base")
#endif
extern doublereal dlapy2_(doublereal *x, doublereal *y);
/*
f2c knows the exact rules for precedence, and so omits parentheses where not
strictly necessary. Since this is generated code, we don't really care if
it's readable, and we know what is written is correct. So don't warn about
them.
*/
#if defined(__GNUC__)
#pragma GCC diagnostic ignored "-Wparentheses"
#endif
'''
class FortranRoutine:
"""Wrapper for a Fortran routine in a file.
"""
type = 'generic'
def __init__(self, name=None, filename=None):
self.filename = filename
if name is None:
root, ext = os.path.splitext(filename)
name = root
self.name = name
self._dependencies = None
def dependencies(self):
if self._dependencies is None:
deps = fortran.getDependencies(self.filename)
self._dependencies = [d.lower() for d in deps]
return self._dependencies
def __repr__(self):
return "FortranRoutine({!r}, filename={!r})".format(self.name,
self.filename)
class UnknownFortranRoutine(FortranRoutine):
"""Wrapper for a Fortran routine for which the corresponding file
is not known.
"""
type = 'unknown'
def __init__(self, name):
FortranRoutine.__init__(self, name=name, filename='<unknown>')
def dependencies(self):
return []
class FortranLibrary:
"""Container for a bunch of Fortran routines.
"""
def __init__(self, src_dirs):
self._src_dirs = src_dirs
self.names_to_routines = {}
def _findRoutine(self, rname):
# 将输入的例程名转换为小写
rname = rname.lower()
# 遍历源目录列表中的每一个目录
for s in self._src_dirs:
# 构建可能的Fortran例程文件名
ffilename = os.path.join(s, rname + '.f')
# 如果文件存在，返回新创建的Fortran例程对象
if os.path.exists(ffilename):
return self._newFortranRoutine(rname, ffilename)
# 如果未找到文件，返回一个未知的Fortran例程对象
return UnknownFortranRoutine(rname)
def _newFortranRoutine(self, rname, filename):
# 创建并返回一个新的Fortran例程对象
return FortranRoutine(rname, filename)
def addIgnorableRoutine(self, rname):
"""Add a routine that we don't want to consider when looking at
dependencies.
"""
# 将例程名转换为小写
rname = rname.lower()
# 创建一个未知的Fortran例程对象
routine = UnknownFortranRoutine(rname)
# 将例程名和对象存储到类属性中
self.names_to_routines[rname] = routine
def addRoutine(self, rname):
"""Add a routine to the library.
"""
# 调用getRoutine方法以确保例程存在于库中
self.getRoutine(rname)
def getRoutine(self, rname):
"""Get a routine from the library. Will add if it's not found.
"""
# 创建一个唯一标识符列表
unique = []
# 将例程名转换为小写
rname = rname.lower()
# 从类属性中获取对应例程名的例程对象，若不存在则返回唯一标识符列表
routine = self.names_to_routines.get(rname, unique)
# 如果获取的是唯一标识符列表，则调用_findRoutine方法寻找并添加例程
if routine is unique:
routine = self._findRoutine(rname)
self.names_to_routines[rname] = routine
# 返回找到或者新添加的例程对象
return routine
def allRoutineNames(self):
"""Return the names of all the routines.
"""
# 返回类属性中所有例程名组成的列表
return list(self.names_to_routines.keys())
def allRoutines(self):
"""Return all the routines.
"""
# 返回类属性中所有例程对象组成的列表
return list(self.names_to_routines.values())
def resolveAllDependencies(self):
"""Try to add routines to the library to satisfy all the dependencies
for each routine in the library.
Returns a set of routine names that have the dependencies unresolved.
"""
# 已处理的例程名集合
done_this = set()
# 上次循环未解决的例程名集合
last_todo = set()
while True:
# 获取所有未处理的例程名集合
todo = set(self.allRoutineNames()) - done_this
# 如果当前未处理集合与上次相同，则退出循环
if todo == last_todo:
break
# 遍历每一个未处理的例程名
for rn in todo:
# 获取该例程名对应的例程对象
r = self.getRoutine(rn)
# 获取该例程对象的所有依赖例程名
deps = r.dependencies()
# 为每一个依赖例程名添加到库中
for d in deps:
self.addRoutine(d)
# 将该例程名标记为已处理
done_this.add(rn)
# 更新上次未处理的例程名集合
last_todo = todo
# 返回未解决依赖的例程名集合
return todo
# LapackLibrary 类，继承自 FortranLibrary 类
class LapackLibrary(FortranLibrary):
# 重写父类方法 _newFortranRoutine
def _newFortranRoutine(self, rname, filename):
# 调用父类方法创建新的 FortranRoutine 对象
routine = FortranLibrary._newFortranRoutine(self, rname, filename)
# 根据文件名和函数名设置 routine 的 type 属性
if 'blas' in filename.lower():
routine.type = 'blas'
elif 'install' in filename.lower():
routine.type = 'config'
elif rname.startswith('z'):
routine.type = 'z_lapack'
elif rname.startswith('c'):
routine.type = 'c_lapack'
elif rname.startswith('s'):
routine.type = 's_lapack'
elif rname.startswith('d'):
routine.type = 'd_lapack'
else:
routine.type = 'lapack'
return routine
# 打印描述和给定类型的所有 routines 名称
def printRoutineNames(desc, routines):
print(desc)
for r in routines:
print('\t%s' % r.name)
# 根据给定的 wrapped_routines、ignores 和 lapack_dir 创建 LapackLibrary 对象
def getLapackRoutines(wrapped_routines, ignores, lapack_dir):
# 确定 BLAS 和 LAPACK 源码目录
blas_src_dir = os.path.join(lapack_dir, 'BLAS', 'SRC')
if not os.path.exists(blas_src_dir):
blas_src_dir = os.path.join(lapack_dir, 'blas', 'src')
lapack_src_dir = os.path.join(lapack_dir, 'SRC')
if not os.path.exists(lapack_src_dir):
lapack_src_dir = os.path.join(lapack_dir, 'src')
install_src_dir = os.path.join(lapack_dir, 'INSTALL')
if not os.path.exists(install_src_dir):
install_src_dir = os.path.join(lapack_dir, 'install')
# 创建 LapackLibrary 对象
library = LapackLibrary([install_src_dir, blas_src_dir, lapack_src_dir])
# 将 ignores 中的 routine 添加到忽略列表中
for r in ignores:
library.addIgnorableRoutine(r)
# 将 wrapped_routines 中的 routine 添加到 library 中
for w in wrapped_routines:
library.addRoutine(w)
# 解析所有依赖关系
library.resolveAllDependencies()
return library
# 从 wrapped_routines_file 中获取 wrapped routines 和 ignores
def getWrappedRoutineNames(wrapped_routines_file):
routines = []
ignores = []
with open(wrapped_routines_file) as fo:
for line in fo:
line = line.strip()
if not line or line.startswith('#'):
continue
if line.startswith('IGNORE:'):
line = line[7:].strip()
ig = line.split()
ignores.extend(ig)
else:
routines.append(line)
return routines, ignores
# 类型集合，包含了可能的 routine 类型
types = {'blas', 'lapack', 'd_lapack', 's_lapack', 'z_lapack', 'c_lapack', 'config'}
# 将 library 中每种类型的 routines 名称和依赖关系写入输出目录中的文件
def dumpRoutineNames(library, output_dir):
for typename in {'unknown'} | types:
routines = library.allRoutinesByType(typename)
filename = os.path.join(output_dir, typename + '_routines.lst')
with open(filename, 'w') as fo:
for r in routines:
deps = r.dependencies()
fo.write('%s: %s\n' % (r.name, ' '.join(deps)))
# 将 routines 中的所有源文件内容合并并写入 output_file
def concatenateRoutines(routines, output_file):
with open(output_file, 'w') as output_fo:
for r in routines:
with open(r.filename) as fo:
source = fo.read()
output_fo.write(source)
class F2CError(Exception):
pass
# 将 Fortran 文件名中的反斜杠替换为正斜杠，以适应不同操作系统路径格式
def runF2C(fortran_filename, output_dir):
fortran_filename = fortran_filename.replace('\\', '/')
try:
# 使用 subprocess 调用 f2c 转换 Fortran 文件为 C 文件，并指定输出目录
subprocess.check_call(
["f2c"] + F2C_ARGS + ['-d', output_dir, fortran_filename]
)
except subprocess.CalledProcessError:
# 如果调用出错，则抛出自定义异常 F2CError
raise F2CError
# 清理指定 C 文件中的源码，通过 clapack_scrub 模块进行处理
def scrubF2CSource(c_file):
with open(c_file) as fo:
source = fo.read()
# 调用 clapack_scrub 模块的 scrubSource 函数清理源码，并输出详细信息
source = clapack_scrub.scrubSource(source, verbose=True)
with open(c_file, 'w') as fo:
# 将清理后的源码写回原文件，并添加标准头部 HEADER
fo.write(HEADER)
fo.write(source)
# 确保指定的可执行文件名在系统路径中可找到，否则抛出 SystemExit 异常
def ensure_executable(name):
try:
which(name)
except Exception:
raise SystemExit(name + ' not found')
# 创建 LAPACK 符号重命名头文件，以避免符号冲突，并加入 BLAS/LAPACK 和 f2c 的重命名
def create_name_header(output_dir):
routine_re = re.compile(r'^ (subroutine|.* function)\s+(\w+)\(.*$',
re.I)
extern_re = re.compile(r'^extern [a-z]+ ([a-z0-9_]+)\(.*$')
# BLAS/LAPACK 符号集合
symbols = set(['xerbla'])
for fn in os.listdir(output_dir):
fn = os.path.join(output_dir, fn)
if not fn.endswith('.f'):
continue
with open(fn) as f:
for line in f:
m = routine_re.match(line)
if m:
# 提取并添加符号到集合中
symbols.add(m.group(2).lower())
# f2c 符号集合
f2c_symbols = set()
with open('f2c.h') as f:
for line in f:
m = extern_re.match(line)
if m:
# 提取并添加 f2c 符号到集合中
f2c_symbols.add(m.group(1))
with open(os.path.join(output_dir, 'lapack_lite_names.h'), 'w') as f:
f.write(HEADER_BLURB)
f.write(
"/*\n"
" * This file renames all BLAS/LAPACK and f2c symbols to avoid\n"
" * dynamic symbol name conflicts, in cases where e.g.\n"
" * integer sizes do not match with 'standard' ABI.\n"
" */\n")
# 重命名 BLAS/LAPACK 符号并写入文件
for name in sorted(symbols):
f.write("#define %s_ BLAS_FUNC(%s)\n" % (name, name))
# 也重命名 f2c 导出的符号并写入文件
f.write("\n"
"/* Symbols exported by f2c.c */\n")
for name in sorted(f2c_symbols):
f.write("#define %s numpy_lapack_lite_%s\n" % (name, name))
# 主程序入口，执行 LAPACK 转换和符号重命名操作
def main():
if len(sys.argv) != 3:
# 如果命令行参数不为3个，打印帮助文档并返回
print(__doc__)
return
# 确保系统能找到 patch 和 f2c 可执行文件
ensure_executable('f2c')
ensure_executable('patch')
# 读取命令行参数
wrapped_routines_file = sys.argv[1]
lapack_src_dir = sys.argv[2]
output_dir = os.path.join(os.path.dirname(__file__), 'build')
# 删除旧的输出目录并创建新的输出目录
shutil.rmtree(output_dir, ignore_errors=True)
os.makedirs(output_dir)
# 获取需要转换的 LAPACK 函数列表和忽略列表
wrapped_routines, ignores = getWrappedRoutineNames(wrapped_routines_file)
# 获取 LAPACK 函数库
library = getLapackRoutines(wrapped_routines, ignores, lapack_src_dir)
# 将 LAPACK 函数信息导出到指定的输出目录
dumpRoutineNames(library, output_dir)
# 遍历类型列表中的每个类型名
for typename in types:
# 构建对应的 Fortran 文件路径
fortran_file = os.path.join(output_dir, 'f2c_%s.f' % typename)
# 根据 Fortran 文件路径生成对应的 C 文件路径
c_file = fortran_file[:-2] + '.c'
# 打印正在创建的 C 文件名
print('creating %s ...' % c_file)
# 获取特定类型的所有例程
routines = library.allRoutinesByType(typename)
# 将所有例程连接起来并写入 Fortran 文件
concatenateRoutines(routines, fortran_file)
# 应用补丁文件
patch_file = os.path.basename(fortran_file) + '.patch'
# 如果补丁文件存在，则使用 subprocess 调用 patch 命令进行打补丁操作
if os.path.exists(patch_file):
subprocess.check_call(['patch', '-u', fortran_file, patch_file])
print("Patched {}".format(fortran_file))
try:
# 尝试运行 f2c 工具转换 Fortran 文件为 C 文件
runF2C(fortran_file, output_dir)
except F2CError:
# 如果转换失败，则打印错误信息
print('f2c failed on %s' % fortran_file)
break
# 清理生成的 C 文件的源代码
scrubF2CSource(c_file)
# 检查并应用 C 文件的补丁
c_patch_file = os.path.basename(c_file) + '.patch'
if os.path.exists(c_patch_file):
subprocess.check_call(['patch', '-u', c_file, c_patch_file])
# 打印空行，用于分隔不同文件处理的输出信息
print()
# 创建名称头文件
create_name_header(output_dir)
# 遍历输出目录中的文件
for fname in os.listdir(output_dir):
# 复制所有以 '.c' 结尾的文件或特定的头文件 'lapack_lite_names.h'
if fname.endswith('.c') or fname == 'lapack_lite_names.h':
print('Copying ' + fname)
# 将文件复制到当前脚本所在目录
shutil.copy(
os.path.join(output_dir, fname),
os.path.abspath(os.path.dirname(__file__)),
)
# 如果当前脚本作为主程序运行
if __name__ == '__main__':
# 调用主函数 main()
main()
.\numpy\numpy\linalg\lapack_lite\python_xerbla.c
/*
定义宏以清除 Py_ssize_t 的定义
包含 Python.h 头文件
包含 numpy 库的通用头文件和 nmpy_cblas.h 头文件
*/
/*
从原始手册页：
--------------------------
XERBLA 是 LAPACK 程序的错误处理程序。
如果输入参数具有无效值，则由 LAPACK 程序调用它。
将打印一条消息并停止执行。
不打印消息和停止执行，而是引发一个带有消息的 ValueError 异常。
参数:
-----------
srname: 错误消息中使用的子程序名称，最多六个字符。
结尾的空格会被跳过。
info: 无效参数的编号。
*/
CBLAS_INT BLAS_FUNC(xerbla)(char *srname, CBLAS_INT *info)
{
static const char format[] = "On entry to %.*s" \
" parameter number %d had an illegal value";
char buf[sizeof(format) + 6 + 4]; /* 6 for name, 4 for param. num. */
int len = 0; /* 子程序名称的长度 */
PyGILState_STATE save;
while( len<6 && srname[len]!='\0' )
len++;
while( len && srname[len-1]==' ' )
len--;
// 保证全局解释器锁（GIL）被获取
save = PyGILState_Ensure();
// 使用给定的格式化字符串生成错误消息
PyOS_snprintf(buf, sizeof(buf), format, len, srname, (int)*info);
// 设置 ValueError 异常并传递错误消息
PyErr_SetString(PyExc_ValueError, buf);
// 释放全局解释器锁（GIL）
PyGILState_Release(save);
// 返回 0，表示出错处理函数的返回值
return 0;
}
.\numpy\numpy\linalg\lapack_litemodule.c
/*
* 此模块由Doug Heisterkamp贡献
* Jim Hugunin进行了修改
* Jeff Whitaker进行了更多修改
*/
/* 定义取消过时 API 的宏 */
#define NPY_NO_DEPRECATED_API NPY_API_VERSION
/* 清除以前的 ssize_t 定义，以使用更安全的定义 */
#define PY_SSIZE_T_CLEAN
#include <Python.h>
/* 包含 NumPy 的数组对象头文件 */
#include "numpy/arrayobject.h"
/* 包含 NumPy 的 CBLAS 头文件 */
#include "npy_cblas.h"
/* 定义一个宏，用于生成 BLAS 函数名 */
#define FNAME(name) BLAS_FUNC(name)
/* 定义 Fortran 中使用的整型 */
typedef CBLAS_INT fortran_int;
#ifdef HAVE_BLAS_ILP64
/* 根据不同的数据类型和平台选择对应的 Fortran 整型格式 */
#if NPY_BITSOF_SHORT == 64
#define FINT_PYFMT "h"
#elif NPY_BITSOF_INT == 64
#define FINT_PYFMT "i"
#elif NPY_BITSOF_LONG == 64
#define FINT_PYFMT "l"
#elif NPY_BITSOF_LONGLONG == 64
#define FINT_PYFMT "L"
#else
/* 如果没有找到兼容的 64 位整型大小，则抛出错误 */
#error No compatible 64-bit integer size. \
Please contact NumPy maintainers and give detailed information about your \
compiler and platform, or dont try to use ILP64 BLAS
#endif
#else
/* 默认使用标准的 32 位整型 */
#define FINT_PYFMT "i"
#endif
/* 定义复数类型结构体 */
typedef struct { float r, i; } f2c_complex;
typedef struct { double r, i; } f2c_doublecomplex;
/* typedef long int (*L_fp)(); */
/* 声明外部链接的 BLAS 和 LAPACK 函数 */
extern fortran_int FNAME(dgelsd)(fortran_int *m, fortran_int *n, fortran_int *nrhs,
double a[], fortran_int *lda, double b[], fortran_int *ldb,
double s[], double *rcond, fortran_int *rank,
double work[], fortran_int *lwork, fortran_int iwork[], fortran_int *info);
extern fortran_int FNAME(zgelsd)(fortran_int *m, fortran_int *n, fortran_int *nrhs,
f2c_doublecomplex a[], fortran_int *lda,
f2c_doublecomplex b[], fortran_int *ldb,
double s[], double *rcond, fortran_int *rank,
f2c_doublecomplex work[], fortran_int *lwork,
double rwork[], fortran_int iwork[], fortran_int *info);
extern fortran_int FNAME(dgeqrf)(fortran_int *m, fortran_int *n, double a[], fortran_int *lda,
double tau[], double work[],
fortran_int *lwork, fortran_int *info);
extern fortran_int FNAME(zgeqrf)(fortran_int *m, fortran_int *n, f2c_doublecomplex a[], fortran_int *lda,
f2c_doublecomplex tau[], f2c_doublecomplex work[],
fortran_int *lwork, fortran_int *info);
extern fortran_int FNAME(dorgqr)(fortran_int *m, fortran_int *n, fortran_int *k, double a[], fortran_int *lda,
double tau[], double work[],
fortran_int *lwork, fortran_int *info);
extern fortran_int FNAME(zungqr)(fortran_int *m, fortran_int *n, fortran_int *k, f2c_doublecomplex a[],
fortran_int *lda, f2c_doublecomplex tau[],
f2c_doublecomplex work[], fortran_int *lwork, fortran_int *info);
extern fortran_int FNAME(xerbla)(char *srname, fortran_int *info);
/* 定义 LapackError 的静态 PyObject 指针 */
static PyObject *LapackError;
/* 定义一个宏，用于简化错误处理 */
#define TRY(E) if (!(E)) return NULL
/* 定义一个静态函数，用于检查 Python 对象 */
static int
check_object(PyObject *ob, int t, char *obname,
char *tname, char *funname)
{
# 检查参数 ob 是否为 NumPy 数组，如果不是则抛出错误信息并返回 0
if (!PyArray_Check(ob)) {
PyErr_Format(LapackError,
"Expected an array for parameter %s in lapack_lite.%s",
obname, funname);
return 0;
}
# 检查参数 ob 是否为 C 连续存储的 NumPy 数组，如果不是则抛出错误信息并返回 0
else if (!PyArray_IS_C_CONTIGUOUS((PyArrayObject *)ob)) {
PyErr_Format(LapackError,
"Parameter %s is not contiguous in lapack_lite.%s",
obname, funname);
return 0;
}
# 检查参数 ob 是否为指定类型 t 的 NumPy 数组，如果不是则抛出错误信息并返回 0
else if (!(PyArray_TYPE((PyArrayObject *)ob) == t)) {
PyErr_Format(LapackError,
"Parameter %s is not of type %s in lapack_lite.%s",
obname, tname, funname);
return 0;
}
# 检查参数 ob 是否具有非本机字节顺序，如果是则抛出错误信息并返回 0
else if (PyArray_ISBYTESWAPPED((PyArrayObject *)ob)) {
PyErr_Format(LapackError,
"Parameter %s has non-native byte order in lapack_lite.%s",
obname, funname);
return 0;
}
# 如果所有条件都满足，则返回 1 表示通过参数检查
else {
return 1;
}
# 定义宏 CHDATA，将指针 p 强制转换为 char* 类型，并获取其指向的 NumPy 数组数据的指针
#define CHDATA(p) ((char *) PyArray_DATA((PyArrayObject *)p))
# 定义宏 SHDATA，将指针 p 强制转换为 short int* 类型，并获取其指向的 NumPy 数组数据的指针
#define SHDATA(p) ((short int *) PyArray_DATA((PyArrayObject *)p))
# 定义宏 DDATA，将指针 p 强制转换为 double* 类型，并获取其指向的 NumPy 数组数据的指针
#define DDATA(p) ((double *) PyArray_DATA((PyArrayObject *)p))
# 定义宏 FDATA，将指针 p 强制转换为 float* 类型，并获取其指向的 NumPy 数组数据的指针
#define FDATA(p) ((float *) PyArray_DATA((PyArrayObject *)p))
# 定义宏 CDATA，将指针 p 强制转换为 f2c_complex* 类型，并获取其指向的 NumPy 数组数据的指针
#define CDATA(p) ((f2c_complex *) PyArray_DATA((PyArrayObject *)p))
# 定义宏 ZDATA，将指针 p 强制转换为 f2c_doublecomplex* 类型，并获取其指向的 NumPy 数组数据的指针
#define ZDATA(p) ((f2c_doublecomplex *) PyArray_DATA((PyArrayObject *)p))
# 定义宏 IDATA，将指针 p 强制转换为 fortran_int* 类型，并获取其指向的 NumPy 数组数据的指针
#define IDATA(p) ((fortran_int *) PyArray_DATA((PyArrayObject *)p))
# 定义 lapack_lite_dgelsd 函数，接收 Python 的调用，执行 dgelsd LAPACK 求解操作
static PyObject *
lapack_lite_dgelsd(PyObject *NPY_UNUSED(self), PyObject *args)
{
fortran_int lapack_lite_status; // LAPACK 操作返回状态
fortran_int m; // 矩阵 A 的行数
fortran_int n; // 矩阵 A 的列数
fortran_int nrhs; // 矩阵 B 的列数
PyObject *a; // NumPy 数组对象，存储矩阵 A 的数据
fortran_int lda; // 矩阵 A 的列数
PyObject *b; // NumPy 数组对象，存储矩阵 B 的数据
fortran_int ldb; // 矩阵 B 的列数
PyObject *s; // NumPy 数组对象，存储奇异值分解结果 S 的数据
double rcond; // 控制奇异值的截断参数
fortran_int rank; // 估计的矩阵 A 的秩
PyObject *work; // NumPy 数组对象，提供给 LAPACK 的工作空间
PyObject *iwork; // NumPy 数组对象，提供给 LAPACK 的整数工作空间
fortran_int lwork; // 工作空间的长度
fortran_int info; // LAPACK 操作的信息
TRY(PyArg_ParseTuple(args,
(FINT_PYFMT FINT_PYFMT FINT_PYFMT "O" FINT_PYFMT "O"
FINT_PYFMT "O" "d" FINT_PYFMT "O" FINT_PYFMT "O"
FINT_PYFMT ":dgelsd"),
&m,&n,&nrhs,&a,&lda,&b,&ldb,&s,&rcond,
&rank,&work,&lwork,&iwork,&info));
TRY(check_object(a,NPY_DOUBLE,"a","NPY_DOUBLE","dgelsd")); // 检查参数 a 是否为 NPY_DOUBLE 类型的 NumPy 数组
TRY(check_object(b,NPY_DOUBLE,"b","NPY_DOUBLE","dgelsd")); // 检查参数 b 是否为 NPY_DOUBLE 类型的 NumPy 数组
TRY(check_object(s,NPY_DOUBLE,"s","NPY_DOUBLE","dgelsd")); // 检查参数 s 是否为 NPY_DOUBLE 类型的 NumPy 数组
TRY(check_object(work,NPY_DOUBLE,"work","NPY_DOUBLE","dgelsd")); // 检查参数 work 是否为 NPY_DOUBLE 类型的 NumPy 数组
#ifndef NPY_UMATH_USE_BLAS64_
TRY(check_object(iwork,NPY_INT,"iwork","NPY_INT","dgelsd")); // 检查参数 iwork 是否为 NPY_INT 类型的 NumPy 数组
#else
TRY(check_object(iwork,NPY_INT64,"iwork","NPY_INT64","dgelsd")); // 检查参数 iwork 是否为 NPY_INT64 类型的 NumPy 数组
#endif
// 调用 LAPACK 中的 dgelsd 函数，进行最小二乘法求解
lapack_lite_status =
FNAME(dgelsd)(&m,&n,&nrhs,DDATA(a),&lda,DDATA(b),&ldb,
DDATA(s),&rcond,&rank,DDATA(work),&lwork,
IDATA(iwork),&info);
if (PyErr_Occurred()) {
return NULL;
}
// 构建返回值，包括 dgelsd 操作的结果信息
return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
",s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
",s:d,s:" FINT_PYFMT ",s:" FINT_PYFMT
",s:" FINT_PYFMT "}"),
"dgelsd_",lapack_lite_status,"m",m,"n",n,"nrhs",nrhs,
"lda",lda,"ldb",ldb,"rcond",rcond,"rank",rank,
"lwork",lwork,"info",info);
}
{
fortran_int lapack_lite_status; // 定义 LAPACK 函数返回状态变量
fortran_int m, n, lwork; // 定义矩阵维度 m, n 和工作数组大小 lwork
PyObject *a, *tau, *work; // 定义 Python 对象指针，分别用于矩阵 A、TAU 和工作数组
fortran_int lda; // 定义矩阵 A 的 leading dimension
fortran_int info; // 定义 LAPACK 函数信息变量
TRY(PyArg_ParseTuple(args,
(FINT_PYFMT FINT_PYFMT "O" FINT_PYFMT "OO"
FINT_PYFMT FINT_PYFMT ":dgeqrf"),
&m, &n, &a, &lda, &tau, &work, &lwork, &info));
// 解析传入的 Python 元组参数，获取 m, n, a, lda, tau, work, lwork 和 info 的值
/* check objects and convert to right storage order */
TRY(check_object(a, NPY_DOUBLE, "a", "NPY_DOUBLE", "dgeqrf"));
TRY(check_object(tau, NPY_DOUBLE, "tau", "NPY_DOUBLE", "dgeqrf"));
TRY(check_object(work, NPY_DOUBLE, "work", "NPY_DOUBLE", "dgeqrf"));
// 检查对象类型并确保存储顺序正确，针对矩阵 a、tau 和 work
lapack_lite_status =
FNAME(dgeqrf)(&m, &n, DDATA(a), &lda, DDATA(tau),
DDATA(work), &lwork, &info);
// 调用 LAPACK 的 dgeqrf 函数进行 QR 分解计算
if (PyErr_Occurred()) {
return NULL;
}
// 检查是否有 Python 异常发生，如有则返回空指针
return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
",s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT "}"),
"dgeqrf_",
lapack_lite_status, "m", m, "n", n, "lda", lda,
"lwork", lwork, "info", info);
// 构建并返回包含函数执行结果的 Python 元组对象
}
static PyObject *
lapack_lite_dorgqr(PyObject *NPY_UNUSED(self), PyObject *args)
{
fortran_int lapack_lite_status;
fortran_int m, n, k, lwork;
PyObject *a, *tau, *work;
fortran_int lda;
fortran_int info;
TRY(PyArg_ParseTuple(args,
(FINT_PYFMT FINT_PYFMT FINT_PYFMT "O"
FINT_PYFMT "OO" FINT_PYFMT FINT_PYFMT
":dorgqr"),
&m, &n, &k, &a, &lda, &tau, &work, &lwork, &info));
// 解析传入的 Python 元组参数，获取 m, n, k, a, lda, tau, work, lwork 和 info 的值
TRY(check_object(a, NPY_DOUBLE, "a", "NPY_DOUBLE", "dorgqr"));
TRY(check_object(tau, NPY_DOUBLE, "tau", "NPY_DOUBLE", "dorgqr"));
TRY(check_object(work, NPY_DOUBLE, "work", "NPY_DOUBLE", "dorgqr"));
// 检查对象类型并确保存储顺序正确，针对矩阵 a、tau 和 work
lapack_lite_status =
FNAME(dorgqr)(&m, &n, &k, DDATA(a), &lda, DDATA(tau), DDATA(work),
&lwork, &info);
// 调用 LAPACK 的 dorgqr 函数进行 QR 分解计算
if (PyErr_Occurred()) {
return NULL;
}
// 检查是否有 Python 异常发生，如有则返回空指针
return Py_BuildValue("{s:i,s:i}", "dorgqr_", lapack_lite_status,
"info", info);
// 构建并返回包含函数执行结果的 Python 字典对象
}
static PyObject *
lapack_lite_zgelsd(PyObject *NPY_UNUSED(self), PyObject *args)
{
fortran_int lapack_lite_status;
fortran_int m;
fortran_int n;
fortran_int nrhs;
PyObject *a;
fortran_int lda;
PyObject *b;
fortran_int ldb;
PyObject *s;
double rcond;
fortran_int rank;
PyObject *work;
fortran_int lwork;
PyObject *rwork;
PyObject *iwork;
fortran_int info;
# 尝试解析传入的参数元组，根据指定格式字符串匹配参数，对应关系如下：
# m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, rwork, iwork, info
TRY(PyArg_ParseTuple(args,
(FINT_PYFMT FINT_PYFMT FINT_PYFMT "O" FINT_PYFMT
"O" FINT_PYFMT "Od" FINT_PYFMT "O" FINT_PYFMT
"OO" FINT_PYFMT ":zgelsd"),
&m,&n,&nrhs,&a,&lda,&b,&ldb,&s,&rcond,
&rank,&work,&lwork,&rwork,&iwork,&info));
# 尝试检查参数对象 a，确保其为 NPY_CDOUBLE 类型，用于 zgelsd 函数
TRY(check_object(a,NPY_CDOUBLE,"a","NPY_CDOUBLE","zgelsd"));
# 尝试检查参数对象 b，确保其为 NPY_CDOUBLE 类型，用于 zgelsd 函数
TRY(check_object(b,NPY_CDOUBLE,"b","NPY_CDOUBLE","zgelsd"));
# 尝试检查参数对象 s，确保其为 NPY_DOUBLE 类型，用于 zgelsd 函数
TRY(check_object(s,NPY_DOUBLE,"s","NPY_DOUBLE","zgelsd"));
# 尝试检查参数对象 work，确保其为 NPY_CDOUBLE 类型，用于 zgelsd 函数
TRY(check_object(work,NPY_CDOUBLE,"work","NPY_CDOUBLE","zgelsd"));
# 尝试检查参数对象 rwork，确保其为 NPY_DOUBLE 类型，用于 zgelsd 函数
TRY(check_object(rwork,NPY_DOUBLE,"rwork","NPY_DOUBLE","zgelsd"));
#ifndef NPY_UMATH_USE_BLAS64_
// 如果未定义 NPY_UMATH_USE_BLAS64_，调用 check_object 函数检查 iwork 对象，期望其类型为 NPY_INT，用于函数 zgelsd
TRY(check_object(iwork,NPY_INT,"iwork","NPY_INT","zgelsd"));
#else
// 如果定义了 NPY_UMATH_USE_BLAS64_，调用 check_object 函数检查 iwork 对象，期望其类型为 NPY_INT64，用于函数 zgelsd
TRY(check_object(iwork,NPY_INT64,"iwork","NPY_INT64","zgelsd"));
#endif
// 调用 LAPACK 函数 zgelsd 执行最小二乘解，对矩阵进行奇异值分解
lapack_lite_status =
FNAME(zgelsd)(&m,&n,&nrhs,ZDATA(a),&lda,ZDATA(b),&ldb,DDATA(s),&rcond,
&rank,ZDATA(work),&lwork,DDATA(rwork),IDATA(iwork),&info);
// 检查是否发生了 Python 异常，若有则返回空指针
if (PyErr_Occurred()) {
return NULL;
}
// 返回 LAPACK 函数 zgelsd 的执行结果和相关参数的 Python 对象表示
return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
",s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
",s:" FINT_PYFMT ",s:" FINT_PYFMT ",s:" FINT_PYFMT
"}"),
"zgelsd_",
lapack_lite_status,"m",m,"n",n,"nrhs",nrhs,"lda",lda,
"ldb",ldb,"rank",rank,"lwork",lwork,"info",info);
}
static PyObject *
lapack_lite_zgeqrf(PyObject *NPY_UNUSED(self), PyObject *args)
{
fortran_int lapack_lite_status;
fortran_int m, n, lwork;
PyObject *a, *tau, *work;
fortran_int lda;
fortran_int info;
// 尝试解析 Python 元组参数 args，解析后的参数依次赋值给 m, n, a, lda, tau, work, lwork, info
TRY(PyArg_ParseTuple(args,
(FINT_PYFMT FINT_PYFMT "O" FINT_PYFMT "OO"
FINT_PYFMT "" FINT_PYFMT ":zgeqrf"),
&m,&n,&a,&lda,&tau,&work,&lwork,&info));
// 检查并转换对象 a, tau, work 的存储顺序为 NPY_CDOUBLE
TRY(check_object(a,NPY_CDOUBLE,"a","NPY_CDOUBLE","zgeqrf"));
TRY(check_object(tau,NPY_CDOUBLE,"tau","NPY_CDOUBLE","zgeqrf"));
TRY(check_object(work,NPY_CDOUBLE,"work","NPY_CDOUBLE","zgeqrf"));
// 调用 LAPACK 函数 zgeqrf 执行 QR 分解
lapack_lite_status =
FNAME(zgeqrf)(&m, &n, ZDATA(a), &lda, ZDATA(tau), ZDATA(work),
&lwork, &info);
// 检查是否发生了 Python 异常，若有则返回空指针
if (PyErr_Occurred()) {
return NULL;
}
// 返回 LAPACK 函数 zgeqrf 的执行结果和相关参数的 Python 对象表示
return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT
",s:" FINT_PYFMT ",s:" FINT_PYFMT
",s:" FINT_PYFMT ",s:" FINT_PYFMT "}"),
"zgeqrf_",lapack_lite_status,"m",m,"n",n,"lda",lda,"lwork",lwork,"info",info);
}
static PyObject *
lapack_lite_zungqr(PyObject *NPY_UNUSED(self), PyObject *args)
{
// lapack_lite 模块中的 zungqr 函数实现
fortran_int lapack_lite_status; // Lapack 函数调用返回状态
fortran_int m, n, k, lwork; // 整型变量 m, n, k, lwork
PyObject *a, *tau, *work; // Python 对象 a, tau, work
fortran_int lda; // Lapack 函数调用中的 lda 参数
fortran_int info; // Lapack 函数调用返回的信息状态
// 解析 Python 函数参数，格式化为 Lapack 函数调用的输入
TRY(PyArg_ParseTuple(args,
(FINT_PYFMT FINT_PYFMT FINT_PYFMT "O"
FINT_PYFMT "OO" FINT_PYFMT "" FINT_PYFMT
":zungqr"),
&m, &n, &k, &a, &lda, &tau, &work, &lwork, &info));
// 检查 Python 对象 a, tau, work 的类型是否符合预期
TRY(check_object(a,NPY_CDOUBLE,"a","NPY_CDOUBLE","zungqr"));
TRY(check_object(tau,NPY_CDOUBLE,"tau","NPY_CDOUBLE","zungqr"));
TRY(check_object(work,NPY_CDOUBLE,"work","NPY_CDOUBLE","zungqr"));
// 调用 Lapack 中的 zungqr 函数进行 QR 分解
lapack_lite_status =
FNAME(zungqr)(&m, &n, &k, ZDATA(a), &lda, ZDATA(tau), ZDATA(work),
&lwork, &info);
// 检查是否有 Python 异常发生，若有则返回空指针
if (PyErr_Occurred()) {
return NULL;
}
// 构建 Python 对象返回结果，包括 zungqr 函数的状态和信息
return Py_BuildValue(("{s:" FINT_PYFMT ",s:" FINT_PYFMT "}"),
"zungqr_",lapack_lite_status,
"info",info);
}
static PyObject *
lapack_lite_xerbla(PyObject *NPY_UNUSED(self), PyObject *args)
{
// 初始化 Lapack 错误信息的状态为 -1
fortran_int info = -1;
NPY_BEGIN_THREADS_DEF; // 定义多线程使用的变量
NPY_BEGIN_THREADS; // 开始多线程执行
// 调用 Lapack 错误处理函数 xerbla
FNAME(xerbla)("test", &info);
NPY_END_THREADS; // 结束多线程执行
// 检查是否有 Python 异常发生，若有则返回空指针
if (PyErr_Occurred()) {
return NULL;
}
// 返回 None
Py_RETURN_NONE;
}
#define STR(x) #x
#define lameth(name) {STR(name), lapack_lite_##name, METH_VARARGS, NULL}
static struct PyMethodDef lapack_lite_module_methods[] = {
lameth(dgelsd), // lapack_lite 模块中的 dgelsd 函数
lameth(dgeqrf), // lapack_lite 模块中的 dgeqrf 函数
lameth(dorgqr), // lapack_lite 模块中的 dorgqr 函数
lameth(zgelsd), // lapack_lite 模块中的 zgelsd 函数
lameth(zgeqrf), // lapack_lite 模块中的 zgeqrf 函数
lameth(zungqr), // lapack_lite 模块中的 zungqr 函数
lameth(xerbla), // lapack_lite 模块中的 xerbla 函数
{ NULL,NULL,0, NULL} // 方法定义结束的标志
};
static struct PyModuleDef moduledef = {
PyModuleDef_HEAD_INIT, // 初始化 Python 模块定义头部
"lapack_lite", // 模块名为 lapack_lite
NULL, // 模块文档字符串为空
-1, // 模块状态为 -1，表示模块全局变量和线程的隔离
lapack_lite_module_methods, // 模块中包含的方法列表
NULL, // 模块的全局变量为 NULL
NULL, // 模块的 slot 常规方法为 NULL
NULL, // 模块的特殊内存分配器为 NULL
NULL // 模块的清理函数为 NULL
};
/* 模块初始化函数 */
PyMODINIT_FUNC PyInit_lapack_lite(void)
{
PyObject *m,*d; // Python 对象 m 和 d
m = PyModule_Create(&moduledef); // 创建名为 lapack_lite 的 Python 模块对象
if (m == NULL) {
return NULL;
}
import_array(); // 导入 NumPy 的数组对象
d = PyModule_GetDict(m); // 获取模块 m 的字典对象
LapackError = PyErr_NewException("lapack_lite.LapackError", NULL, NULL);
PyDict_SetItemString(d, "LapackError", LapackError); // 将 LapackError 异常对象加入模块字典
#ifdef HAVE_BLAS_ILP64
PyDict_SetItemString(d, "_ilp64", Py_True); // 设置模块字典中的 _ilp64 为 Py_True
#else
PyDict_SetItemString(d, "_ilp64", Py_False); // 设置模块字典中的 _ilp64 为 Py_False
#endif
return m; // 返回创建的 Python 模块对象
}
.\numpy\numpy\linalg\linalg.py
# 定义一个特殊的方法 __getattr__，用于在获取不存在的属性时执行
def __getattr__(attr_name):
# 导入警告模块，用于输出警告信息
import warnings
# 从 numpy.linalg 模块中导入 _linalg 对象
from numpy.linalg import _linalg
# 尝试获取 _linalg 对象的属性 attr_name
ret = getattr(_linalg, attr_name, None)
# 如果未找到该属性，则引发 AttributeError 异常
if ret is None:
raise AttributeError(
f"module 'numpy.linalg.linalg' has no attribute {attr_name}")
# 发出警告，说明 numpy.linalg.linalg 已被标记为私有，推荐使用 numpy.linalg._linalg 替代
warnings.warn(
"The numpy.linalg.linalg has been made private and renamed to "
"numpy.linalg._linalg. All public functions exported by it are "
f"available from numpy.linalg. Please use numpy.linalg.{attr_name} "
"instead.",
DeprecationWarning,
stacklevel=3
)
# 返回获取到的属性对象
return ret
.\numpy\numpy\linalg\tests\test_deprecations.py
"""Test deprecation and future warnings.
"""
# 导入 numpy 库，简写为 np
import numpy as np
# 从 numpy.testing 模块中导入 assert_warns 函数
from numpy.testing import assert_warns
# 定义测试函数 test_qr_mode_full_future_warning
def test_qr_mode_full_future_warning():
"""Check mode='full' FutureWarning.
In numpy 1.8 the mode options 'full' and 'economic' in linalg.qr were
deprecated. The release date will probably be sometime in the summer
of 2013.
"""
# 创建一个 2x2 的单位矩阵 a
a = np.eye(2)
# 断言调用 np.linalg.qr 函数时会产生 DeprecationWarning，使用 mode='full'
assert_warns(DeprecationWarning, np.linalg.qr, a, mode='full')
# 断言调用 np.linalg.qr 函数时会产生 DeprecationWarning，使用 mode='f'
assert_warns(DeprecationWarning, np.linalg.qr, a, mode='f')
# 断言调用 np.linalg.qr 函数时会产生 DeprecationWarning，使用 mode='economic'
assert_warns(DeprecationWarning, np.linalg.qr, a, mode='economic')
# 断言调用 np.linalg.qr 函数时会产生 DeprecationWarning，使用 mode='e'
assert_warns(DeprecationWarning, np.linalg.qr, a, mode='e')
.\numpy\numpy\linalg\tests\test_linalg.py
# 导入所需的标准库和第三方库模块
import os
import sys
import itertools
import traceback
import textwrap
import subprocess
import pytest
# 导入 NumPy 库及其部分子模块和函数
import numpy as np
from numpy import array, single, double, csingle, cdouble, dot, identity, matmul
from numpy._core import swapaxes
from numpy.exceptions import AxisError
from numpy import multiply, atleast_2d, inf, asarray
from numpy import linalg
from numpy.linalg import matrix_power, norm, matrix_rank, multi_dot, LinAlgError
from numpy.linalg._linalg import _multi_dot_matrix_chain_order
from numpy.testing import (
assert_, assert_equal, assert_raises, assert_array_equal,
assert_almost_equal, assert_allclose, suppress_warnings,
assert_raises_regex, HAS_LAPACK64, IS_WASM
)
try:
import numpy.linalg.lapack_lite
except ImportError:
# 在没有 BLAS/LAPACK 的情况下可能会导致 ImportError，需要处理这种情况
# 应确保整个测试套件不会因此而中断
pass
def consistent_subclass(out, in_):
# 对于 ndarray 的子类输入，输出应保持相同的子类类型
# （非 ndarray 输入将转换为 ndarray）。
return type(out) is (type(in_) if isinstance(in_, np.ndarray)
else np.ndarray)
old_assert_almost_equal = assert_almost_equal
def assert_almost_equal(a, b, single_decimal=6, double_decimal=12, **kw):
# 根据输入数组的数据类型确定精度，执行数值近似相等的断言
if asarray(a).dtype.type in (single, csingle):
decimal = single_decimal
else:
decimal = double_decimal
old_assert_almost_equal(a, b, decimal=decimal, **kw)
def get_real_dtype(dtype):
# 根据复数数据类型确定其对应的实数数据类型
return {single: single, double: double,
csingle: single, cdouble: double}[dtype]
def get_complex_dtype(dtype):
# 根据实数数据类型确定其对应的复数数据类型
return {single: csingle, double: cdouble,
csingle: csingle, cdouble: cdouble}[dtype]
def get_rtol(dtype):
# 根据数据类型选择一个安全的相对误差阈值
if dtype in (single, csingle):
return 1e-5
else:
return 1e-11
# 用于对测试进行分类的标签集合
all_tags = {
'square', 'nonsquare', 'hermitian', # 互斥的基本类别
'generalized', 'size-0', 'strided' # 可选的附加类别
}
class LinalgCase:
def __init__(self, name, a, b, tags=set()):
"""
用于测试用例的一组参数，包括标识名称、操作数 a 和 b，以及一组标签用于过滤测试
"""
assert_(isinstance(name, str))
self.name = name
self.a = a
self.b = b
self.tags = frozenset(tags) # 使用 frozenset 防止标签共享
def check(self, do):
"""
在此测试用例上运行函数 `do`，并扩展参数
"""
do(self.a, self.b, tags=self.tags)
def __repr__(self):
return f'<LinalgCase: {self.name}>'
def apply_tag(tag, cases):
"""
将给定的标签（字符串）添加到每个测试用例（LinalgCase 对象列表）中
"""
assert tag in all_tags, "Invalid tag"
# 遍历列表中的每个案例对象
for case in cases:
# 将当前案例对象的标签与给定的标签合并
case.tags = case.tags | {tag}
# 返回更新后的案例对象列表
return cases
#
# Base test cases
#
# 设置随机种子，以便重现结果
np.random.seed(1234)
# 创建一个空列表来存储测试用例
CASES = []
# square test cases
# 向 CASES 列表添加“square”标签的测试用例
CASES += apply_tag('square', [
# 创建一个 LinalgCase 对象，包含单精度数组和预期结果
LinalgCase("single",
array([[1., 2.], [3., 4.]], dtype=single),
array([2., 1.], dtype=single)),
LinalgCase("double",
array([[1., 2.], [3., 4.]], dtype=double),
array([2., 1.], dtype=double)),
LinalgCase("double_2",
array([[1., 2.], [3., 4.]], dtype=double),
array([[2., 1., 4.], [3., 4., 6.]], dtype=double)),
LinalgCase("csingle",
array([[1. + 2j, 2 + 3j], [3 + 4j, 4 + 5j]], dtype=csingle),
array([2. + 1j, 1. + 2j], dtype=csingle)),
LinalgCase("cdouble",
array([[1. + 2j, 2 + 3j], [3 + 4j, 4 + 5j]], dtype=cdouble),
array([2. + 1j, 1. + 2j], dtype=cdouble)),
LinalgCase("cdouble_2",
array([[1. + 2j, 2 + 3j], [3 + 4j, 4 + 5j]], dtype=cdouble),
array([[2. + 1j, 1. + 2j, 1 + 3j], [1 - 2j, 1 - 3j, 1 - 6j]], dtype=cdouble)),
LinalgCase("0x0",
np.empty((0, 0), dtype=double),
np.empty((0,), dtype=double),
tags={'size-0'}),
LinalgCase("8x8",
np.random.rand(8, 8),
np.random.rand(8)),
LinalgCase("1x1",
np.random.rand(1, 1),
np.random.rand(1)),
LinalgCase("nonarray",
[[1, 2], [3, 4]],
[2, 1]),
])
# non-square test-cases
# 向 CASES 列表添加“nonsquare”标签的测试用例
CASES += apply_tag('nonsquare', [
LinalgCase("single_nsq_1",
array([[1., 2., 3.], [3., 4., 6.]], dtype=single),
array([2., 1.], dtype=single)),
LinalgCase("single_nsq_2",
array([[1., 2.], [3., 4.], [5., 6.]], dtype=single),
array([2., 1., 3.], dtype=single)),
LinalgCase("double_nsq_1",
array([[1., 2., 3.], [3., 4., 6.]], dtype=double),
array([2., 1.], dtype=double)),
LinalgCase("double_nsq_2",
array([[1., 2.], [3., 4.], [5., 6.]], dtype=double),
array([2., 1., 3.], dtype=double)),
LinalgCase("csingle_nsq_1",
array(
[[1. + 1j, 2. + 2j, 3. - 3j], [3. - 5j, 4. + 9j, 6. + 2j]], dtype=csingle),
array([2. + 1j, 1. + 2j], dtype=csingle)),
LinalgCase("csingle_nsq_2",
array(
[[1. + 1j, 2. + 2j], [3. - 3j, 4. - 9j], [5. - 4j, 6. + 8j]], dtype=csingle),
array([2. + 1j, 1. + 2j, 3. - 3j], dtype=csingle)),
LinalgCase("cdouble_nsq_1",
array(
[[1. + 1j, 2. + 2j, 3. - 3j], [3. - 5j, 4. + 9j, 6. + 2j]], dtype=cdouble),
array([2. + 1j, 1. + 2j], dtype=cdouble)),
LinalgCase("cdouble_nsq_2",
array(
[[1. + 1j, 2. + 2j], [3. - 3j, 4. - 9j], [5. - 4j, 6. + 8j]], dtype=cdouble),
array([2. + 1j, 1. + 2j, 3. - 3j], dtype=cdouble)),
# 创建一个 LinalgCase 对象，用于表示线性代数操作的测试用例，名称为 "cdouble_nsq_1_2"
LinalgCase("cdouble_nsq_1_2",
# 创建一个复数数组，包含两行三列的复数值矩阵
array([[1. + 1j, 2. + 2j, 3. - 3j], [3. - 5j, 4. + 9j, 6. + 2j]], dtype=cdouble),
# 创建一个复数数组，包含两行两列的复数值矩阵
array([[2. + 1j, 1. + 2j], [1 - 1j, 2 - 2j]], dtype=cdouble)),
# 创建一个 LinalgCase 对象，名称为 "cdouble_nsq_2_2"
LinalgCase("cdouble_nsq_2_2",
# 创建一个复数数组，包含三行两列的复数值矩阵
array([[1. + 1j, 2. + 2j], [3. - 3j, 4. - 9j], [5. - 4j, 6. + 8j]], dtype=cdouble),
# 创建一个复数数组，包含三行两列的复数值矩阵
array([[2. + 1j, 1. + 2j], [1 - 1j, 2 - 2j], [1 - 1j, 2 - 2j]], dtype=cdouble)),
# 创建一个 LinalgCase 对象，名称为 "8x11"
LinalgCase("8x11",
# 创建一个 8x11 的随机浮点数数组
np.random.rand(8, 11),
# 创建一个包含 8 个随机浮点数的一维数组
np.random.rand(8)),
# 创建一个 LinalgCase 对象，名称为 "1x5"
LinalgCase("1x5",
# 创建一个 1x5 的随机浮点数数组
np.random.rand(1, 5),
# 创建一个包含 1 个随机浮点数的一维数组
np.random.rand(1)),
# 创建一个 LinalgCase 对象，名称为 "5x1"
LinalgCase("5x1",
# 创建一个 5x1 的随机浮点数数组
np.random.rand(5, 1),
# 创建一个包含 5 个随机浮点数的一维数组
np.random.rand(5)),
# 创建一个 LinalgCase 对象，名称为 "0x4"，带有额外的标签 {'size-0'}
LinalgCase("0x4",
# 创建一个 0x4 的随机浮点数数组
np.random.rand(0, 4),
# 创建一个空的一维数组
np.random.rand(0),
tags={'size-0'}),
# 创建一个 LinalgCase 对象，名称为 "4x0"，带有额外的标签 {'size-0'}
LinalgCase("4x0",
# 创建一个 4x0 的随机浮点数数组
np.random.rand(4, 0),
# 创建一个包含 4 个随机浮点数的一维数组
np.random.rand(4),
tags={'size-0'}),
# 添加测试用例到 hermitian 类别
CASES += apply_tag('hermitian', [
LinalgCase("hsingle",
array([[1., 2.], [2., 1.]], dtype=single),
None),
LinalgCase("hdouble",
array([[1., 2.], [2., 1.]], dtype=double),
None),
LinalgCase("hcsingle",
array([[1., 2 + 3j], [2 - 3j, 1]], dtype=csingle),
None),
LinalgCase("hcdouble",
array([[1., 2 + 3j], [2 - 3j, 1]], dtype=cdouble),
None),
LinalgCase("hempty",
np.empty((0, 0), dtype=double),
None,
tags={'size-0'}),
LinalgCase("hnonarray",
[[1, 2], [2, 1]],
None),
LinalgCase("matrix_b_only",
array([[1., 2.], [2., 1.]]),
None),
LinalgCase("hmatrix_1x1",
np.random.rand(1, 1),
None),
])
# 生成通用测试案例
def _make_generalized_cases():
new_cases = []
for case in CASES:
if not isinstance(case.a, np.ndarray): # 如果 a 不是 ndarray 类型则跳过
continue
a = np.array([case.a, 2 * case.a, 3 * case.a]) # 构造 a 的新数组
if case.b is None: # 如果 b 为空，则设为 None
b = None
elif case.b.ndim == 1:
b = case.b
else:
b = np.array([case.b, 7 * case.b, 6 * case.b])
new_case = LinalgCase(case.name + "_tile3", a, b, tags=case.tags | {'generalized'}) # 创建新的 LinalgCase
new_cases.append(new_case) # 将新的测试案例添加到列表中
a = np.array([case.a] * 2 * 3).reshape((3, 2) + case.a.shape) # 构造 a 的新数组
if case.b is None: # 如果 b 为空，则设为 None
b = None
elif case.b.ndim == 1:
b = np.array([case.b] * 2 * 3 * a.shape[-1]).reshape((3, 2) + case.a.shape[-2:])
else:
b = np.array([case.b] * 2 * 3).reshape((3, 2) + case.b.shape)
new_case = LinalgCase(case.name + "_tile213", a, b, tags=case.tags | {'generalized'}) # 创建新的 LinalgCase
new_cases.append(new_case) # 将新的测试案例添加到列表中
return new_cases # 返回新的测试案例列表
CASES += _make_generalized_cases() # 将通用测试案例添加到 CASES 列表中
def _stride_comb_iter(x):
"""
生成所有轴的步幅的笛卡尔积
"""
if not isinstance(x, np.ndarray): # 如果 x 不是 ndarray 类型
yield x, "nop" # 返回 x 和 "nop"
return
# 初始化步幅组合
stride_set = [(1,)] * x.ndim
stride_set[-1] = (1, 3, -4)
if x.ndim > 1:
stride_set[-2] = (1, 3, -4)
if x.ndim > 2:
stride_set[-3] = (1, -4)
# 使用 itertools.product 生成多个迭代器元组，每个迭代器从 stride_set 中取值
for repeats in itertools.product(*tuple(stride_set)):
# 根据 repeats 计算新的数组形状
new_shape = [abs(a * b) for a, b in zip(x.shape, repeats)]
# 根据 repeats 创建切片对象 slices
slices = tuple([slice(None, None, repeat) for repeat in repeats])
# 创建一个新的数组 xi，具有不同的步幅，但是与 x 具有相同的数据
xi = np.empty(new_shape, dtype=x.dtype)
# 将 xi 视为无符号 32 位整数数组，并填充固定的值 0xdeadbeef
xi.view(np.uint32).fill(0xdeadbeef)
# 根据 slices 切片选择 xi 的数据
xi = xi[slices]
# 将 x 的数据复制到 xi
xi[...] = x
# 将 xi 视为与 x 相同的类别
xi = xi.view(x.__class__)
# 断言 xi 与 x 全部元素相等
assert_(np.all(xi == x))
# 生成一个生成器，返回 xi 和其对应的重复步幅命名
yield xi, "stride_" + "_".join(["%+d" % j for j in repeats])
# 如果可能，生成零步幅的数组
if x.ndim >= 1 and x.shape[-1] == 1:
# 复制 x 的步幅，并将最后一个步幅设为 0
s = list(x.strides)
s[-1] = 0
# 使用 np.lib.stride_tricks.as_strided 创建零步幅的数组 xi
xi = np.lib.stride_tricks.as_strided(x, strides=s)
yield xi, "stride_xxx_0"
if x.ndim >= 2 and x.shape[-2] == 1:
# 复制 x 的步幅，并将倒数第二个步幅设为 0
s = list(x.strides)
s[-2] = 0
xi = np.lib.stride_tricks.as_strided(x, strides=s)
yield xi, "stride_xxx_0_x"
if x.ndim >= 2 and x.shape[:-2] == (1, 1):
# 复制 x 的步幅，并将倒数第一和倒数第二个步幅设为 0
s = list(x.strides)
s[-1] = 0
s[-2] = 0
xi = np.lib.stride_tricks.as_strided(x, strides=s)
yield xi, "stride_xxx_0_0"
def _make_strided_cases():
"""
构造所有可能的步进测试用例，并添加到新的测试用例列表中
"""
new_cases = []
for case in CASES:
for a, a_label in _stride_comb_iter(case.a):
for b, b_label in _stride_comb_iter(case.b):
# 创建带有步进标签的新测试用例对象，并加入到新测试用例列表中
new_case = LinalgCase(case.name + "_" + a_label + "_" + b_label, a, b,
tags=case.tags | {'strided'})
new_cases.append(new_case)
return new_cases
CASES += _make_strided_cases()
"""
将生成的步进测试用例添加到全局测试用例集合中
"""
#
# Test different routines against the above cases
#
class LinalgTestCase:
TEST_CASES = CASES
def check_cases(self, require=set(), exclude=set()):
"""
对每个测试用例执行检查函数，根据require和exclude参数过滤测试用例
"""
for case in self.TEST_CASES:
# 根据require和exclude标签过滤测试用例
if case.tags & require != require:
continue
if case.tags & exclude:
continue
try:
# 执行测试用例的检查函数
case.check(self.do)
except Exception as e:
# 在发生异常时生成详细的错误消息
msg = f'In test case: {case!r}\n\n'
msg += traceback.format_exc()
raise AssertionError(msg) from e
class LinalgSquareTestCase(LinalgTestCase):
def test_sq_cases(self):
"""
运行方阵测试用例，排除广义和零尺寸的情况
"""
self.check_cases(require={'square'},
exclude={'generalized', 'size-0'})
def test_empty_sq_cases(self):
"""
运行空方阵测试用例，排除广义情况但包括零尺寸的情况
"""
self.check_cases(require={'square', 'size-0'},
exclude={'generalized'})
class LinalgNonsquareTestCase(LinalgTestCase):
def test_nonsq_cases(self):
"""
运行非方阵测试用例，排除广义和零尺寸的情况
"""
self.check_cases(require={'nonsquare'},
exclude={'generalized', 'size-0'})
def test_empty_nonsq_cases(self):
"""
运行空非方阵测试用例，排除广义情况但包括零尺寸的情况
"""
self.check_cases(require={'nonsquare', 'size-0'},
exclude={'generalized'})
class HermitianTestCase(LinalgTestCase):
def test_herm_cases(self):
"""
运行厄米特测试用例，排除广义和零尺寸的情况
"""
self.check_cases(require={'hermitian'},
exclude={'generalized', 'size-0'})
def test_empty_herm_cases(self):
"""
运行空厄米特测试用例，排除广义情况但包括零尺寸的情况
"""
self.check_cases(require={'hermitian', 'size-0'},
exclude={'generalized'})
class LinalgGeneralizedSquareTestCase(LinalgTestCase):
@pytest.mark.slow
def test_generalized_sq_cases(self):
"""
运行广义方阵测试用例，包括广义和方阵，但排除零尺寸的情况
"""
self.check_cases(require={'generalized', 'square'},
exclude={'size-0'})
@pytest.mark.slow
def test_generalized_empty_sq_cases(self):
"""
运行空广义方阵测试用例，包括广义、方阵和零尺寸的情况
"""
self.check_cases(require={'generalized', 'square', 'size-0'})
class LinalgGeneralizedNonsquareTestCase(LinalgTestCase):
@pytest.mark.slow
def test_generalized_nonsq_cases(self):
"""
运行广义非方阵测试用例，包括广义和非方阵，但排除零尺寸的情况
"""
self.check_cases(require={'generalized', 'nonsquare'},
exclude={'size-0'})
@pytest.mark.slow
def test_generalized_empty_nonsq_cases(self):
"""
运行空广义非方阵测试用例，包括广义、非方阵和零尺寸的情况
"""
self.check_cases(require={'generalized', 'nonsquare', 'size-0'})
class HermitianGeneralizedTestCase(LinalgTestCase):
@pytest.mark.slow
# 定义一个测试方法，用于测试通用广义厄米特矩阵的情况
def test_generalized_herm_cases(self):
# 调用自定义的方法 check_cases，传入需要的条件集合和排除的条件集合作为参数
self.check_cases(require={'generalized', 'hermitian'},
exclude={'size-0'})
# 使用 pytest 的装饰器标记，将该测试标记为慢速测试
@pytest.mark.slow
# 定义另一个测试方法，用于测试空的通用广义厄米特矩阵的情况
def test_generalized_empty_herm_cases(self):
# 调用自定义的方法 check_cases，传入需要的条件集合、包括空大小的集合和排除的条件集合作为参数
self.check_cases(require={'generalized', 'hermitian', 'size-0'},
exclude={'none'})
# 定义一个函数 `identity_like_generalized`，用于根据输入数组 `a` 的维度返回一个类似单位矩阵的数组
def identity_like_generalized(a):
# 将输入数组转换为 ndarray 类型
a = asarray(a)
# 如果输入数组的维度大于等于3
if a.ndim >= 3:
# 创建一个与输入数组形状相同的空数组，数据类型与输入数组相同
r = np.empty(a.shape, dtype=a.dtype)
# 使用 `identity` 函数填充数组，该函数返回形状为最后两个维度大小的单位矩阵
r[...] = identity(a.shape[-2])
return r
else:
# 如果输入数组的维度小于3，返回一个以输入数组第一个维度大小为形状的单位矩阵
return identity(a.shape[0])
# 定义一个类 `SolveCases`，继承自 `LinalgSquareTestCase` 和 `LinalgGeneralizedSquareTestCase`，用于解决线性代数方程组的测试用例
class SolveCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
# `do` 方法，用于执行测试
# 参数：a - 系数矩阵，b - 常数向量，tags - 标签
def do(self, a, b, tags):
# 使用 `linalg.solve` 求解线性方程组 `a x = b`，返回解向量 `x`
x = linalg.solve(a, b)
# 如果常数向量 `b` 的维度为1
if np.array(b).ndim == 1:
# 计算 `a x`，其中 `b` 视为列向量，结果是 `a x` 的第一维度
adotx = matmul(a, x[..., None])[..., 0]
# 断言 `a x` 扩展到 `b` 的形状后与 `adotx` 相等
assert_almost_equal(np.broadcast_to(b, adotx.shape), adotx)
else:
# 计算 `a x`，直接与 `b` 比较
adotx = matmul(a, x)
# 断言 `a x` 与 `b` 相等
assert_almost_equal(b, adotx)
# 断言 `x` 和 `b` 具有相同的子类
assert_(consistent_subclass(x, b))
# 定义一个类 `TestSolve`，继承自 `SolveCases`，用于测试 `linalg.solve` 的不同数据类型
class TestSolve(SolveCases):
# 测试不同数据类型的参数 `dtype`
@pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
def test_types(self, dtype):
# 创建一个数据类型为 `dtype` 的二维数组 `x`
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
# 断言解 `linalg.solve(x, x)` 的数据类型与输入 `dtype` 相同
assert_equal(linalg.solve(x, x).dtype, dtype)
# 测试输入数组 `a` 和 `b` 的一维情况
def test_1_d(self):
# 定义一个类 `ArraySubclass`，继承自 `np.ndarray`
class ArraySubclass(np.ndarray):
pass
# 创建一个三维数组 `a`
a = np.arange(8).reshape(2, 2, 2)
# 创建一个一维数组 `b`，视图为 `ArraySubclass`
b = np.arange(2).view(ArraySubclass)
# 使用 `linalg.solve` 解方程组 `a x = b`，返回解 `result`
result = linalg.solve(a, b)
# 断言解 `result` 的形状为 (2, 2)
assert result.shape == (2, 2)
# 创建一个二维数组 `b`，视图为 `ArraySubclass`
b = np.arange(4).reshape(2, 2).view(ArraySubclass)
# 使用 `linalg.solve` 解方程组 `a x = b`，返回解 `result`
result = linalg.solve(a, b)
# 断言解 `result` 的形状为 (2, 2, 2)
assert result.shape == (2, 2, 2)
# 创建一个二维数组 `b`，视图为 `ArraySubclass`，形状为 (1, 2)
b = np.arange(2).reshape(1, 2).view(ArraySubclass)
# 断言 `linalg.solve` 抛出 `ValueError` 异常
assert_raises(ValueError, linalg.solve, a, b)
def test_0_size(self):
class ArraySubclass(np.ndarray):
pass
# 创建一个名为 ArraySubclass 的类，继承自 np.ndarray
a = np.arange(8).reshape(2, 2, 2)
# 创建一个 2x2x2 的 NumPy 数组 a
b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)
# 创建一个 1x2x3 的 NumPy 数组 b，并将其视图类型设置为 ArraySubclass
expected = linalg.solve(a, b)[:, 0:0, :]
# 使用 linalg.solve 求解线性方程组 a * x = b，返回结果的第一维在全部情况下、第二维为 0、第三维在全部情况下的部分
result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, :])
# 使用 linalg.solve 求解线性方程组 a[:, 0:0, 0:0] * x = b[:, 0:0, :]，即对 0x0 子数组求解线性方程组
assert_array_equal(result, expected)
# 断言 result 和 expected 数组相等
assert_(isinstance(result, ArraySubclass))
# 断言 result 是 ArraySubclass 类的实例
# 测试非方阵和仅 b 维度为 0 时的错误
assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
# 断言在求解时出现 linalg.LinAlgError 异常，因为 a[:, 0:0, 0:1] 不是方阵
assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])
# 断言在求解时出现 ValueError 异常，因为 b[:, 0:0, :] 维度为 0
# 测试广播错误
b = np.arange(6).reshape(1, 3, 2) # 广播错误
assert_raises(ValueError, linalg.solve, a, b)
# 断言在求解时出现 ValueError 异常，因为 a 和 b 无法广播
assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
# 断言在求解时出现 ValueError 异常，因为 a[0:0] 和 b[0:0] 维度为 0
# 测试 0x0 矩阵的零“单方程”
b = np.arange(2).view(ArraySubclass)
# 创建一个长度为 2 的 ArraySubclass 类型的数组 b
expected = linalg.solve(a, b)[:, 0:0]
# 使用 linalg.solve 求解线性方程组 a * x = b，返回结果的第一维在全部情况下、第二维为 0 的部分
result = linalg.solve(a[:, 0:0, 0:0], b[0:0])
# 使用 linalg.solve 求解线性方程组 a[:, 0:0, 0:0] * x = b[0:0]，即对 0x0 子数组求解线性方程组
assert_array_equal(result, expected)
# 断言 result 和 expected 数组相等
assert_(isinstance(result, ArraySubclass))
# 断言 result 是 ArraySubclass 类的实例
b = np.arange(3).reshape(1, 3)
# 创建一个 1x3 的 NumPy 数组 b
assert_raises(ValueError, linalg.solve, a, b)
# 断言在求解时出现 ValueError 异常，因为 a 和 b 维度不匹配
assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
# 断言在求解时出现 ValueError 异常，因为 a[0:0] 和 b[0:0] 维度为 0
assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b)
# 断言在求解时出现 ValueError 异常，因为 a[:, 0:0, 0:0] 和 b 维度不匹配
class InvCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
# 定义测试用例类 InvCases，继承自 LinalgSquareTestCase 和 LinalgGeneralizedSquareTestCase
def do(self, a, b, tags):
# 实现测试方法 do，接受参数 a, b, tags
a_inv = linalg.inv(a)
# 计算矩阵 a 的逆
assert_almost_equal(matmul(a, a_inv),
identity_like_generalized(a))
# 断言验证 matmul(a, a_inv) 几乎等于 generalized 版本的单位矩阵
assert_(consistent_subclass(a_inv, a))
# 断言验证 a_inv 的类型与 a 一致
class TestInv(InvCases):
# 定义测试类 TestInv，继承自 InvCases
@pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
def test_types(self, dtype):
# 定义测试方法 test_types，参数 dtype 通过 pytest 参数化传入
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
# 创建 dtype 类型的数组 x
assert_equal(linalg.inv(x).dtype, dtype)
# 断言验证 linalg.inv(x) 的数据类型为 dtype
def test_0_size(self):
# 定义测试方法 test_0_size，测试处理各种大小为 0 的数组情况
# 检查所有类型的大小为 0 的数组是否正常工作
class ArraySubclass(np.ndarray):
pass
a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
# 创建大小为 (0, 1, 1) 的 int 类型零数组，视图为 ArraySubclass
res = linalg.inv(a)
# 计算 a 的逆
assert_(res.dtype.type is np.float64)
# 断言验证结果 res 的数据类型是 np.float64
assert_equal(a.shape, res.shape)
# 断言验证 a 和 res 的形状相同
assert_(isinstance(res, ArraySubclass))
# 断言验证 res 是 ArraySubclass 的实例
a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
# 创建大小为 (0, 0) 的 complex64 类型零数组，视图为 ArraySubclass
res = linalg.inv(a)
# 计算 a 的逆
assert_(res.dtype.type is np.complex64)
# 断言验证结果 res 的数据类型是 np.complex64
assert_equal((0,), res.shape)
# 断言验证 res 的形状为 (0,)
assert_(isinstance(res, ArraySubclass))
# 断言验证 res 是 ArraySubclass 的实例
class EigvalsCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
# 定义测试用例类 EigvalsCases，继承自 LinalgSquareTestCase 和 LinalgGeneralizedSquareTestCase
def do(self, a, b, tags):
# 实现测试方法 do，接受参数 a, b, tags
ev = linalg.eigvals(a)
# 计算矩阵 a 的特征值
evalues, evectors = linalg.eig(a)
# 计算矩阵 a 的特征值和特征向量
assert_almost_equal(ev, evalues)
# 断言验证计算得到的特征值 ev 与 linalg.eig 返回的特征值 evalues 几乎相等
class TestEigvals(EigvalsCases):
# 定义测试类 TestEigvals，继承自 EigvalsCases
@pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
def test_types(self, dtype):
# 定义测试方法 test_types，参数 dtype 通过 pytest 参数化传入
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
# 创建 dtype 类型的数组 x
assert_equal(linalg.eigvals(x).dtype, dtype)
# 断言验证 linalg.eigvals(x) 的数据类型为 dtype
x = np.array([[1, 0.5], [-1, 1]], dtype=dtype)
# 创建 dtype 类型的数组 x
assert_equal(linalg.eigvals(x).dtype, get_complex_dtype(dtype))
# 断言验证 linalg.eigvals(x) 的数据类型为复数 dtype 的实部类型
def test_0_size(self):
# 定义测试方法 test_0_size，测试处理各种大小为 0 的数组情况
# 检查所有类型的大小为 0 的数组是否正常工作
class ArraySubclass(np.ndarray):
pass
a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
# 创建大小为 (0, 1, 1) 的 int 类型零数组，视图为 ArraySubclass
res = linalg.eigvals(a)
# 计算 a 的特征值
assert_(res.dtype.type is np.float64)
# 断言验证结果 res 的数据类型是 np.float64
assert_equal((0, 1), res.shape)
# 断言验证 res 的形状为 (0, 1)
# This is just for documentation, it might make sense to change:
assert_(isinstance(res, np.ndarray))
# 断言验证 res 是 np.ndarray 的实例
a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
# 创建大小为 (0, 0) 的 complex64 类型零数组，视图为 ArraySubclass
res = linalg.eigvals(a)
# 计算 a 的特征值
assert_(res.dtype.type is np.complex64)
# 断言验证结果 res 的数据类型是 np.complex64
assert_equal((0,), res.shape)
# 断言验证 res 的形状为 (0,)
# This is just for documentation, it might make sense to change:
assert_(isinstance(res, np.ndarray))
# 断言验证 res 是 np.ndarray 的实例
class EigCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
# 定义测试用例类 EigCases，继承自 LinalgSquareTestCase 和 LinalgGeneralizedSquareTestCase
def do(self, a, b, tags):
# 实现测试方法 do，接受参数 a, b, tags
res = linalg.eig(a)
# 计算矩阵 a 的特征值和特征向量
eigenvalues, eigenvectors = res.eigenvalues, res.eigenvectors
# 获取计算得到的特征值和特征向量
assert_allclose(matmul(a, eigenvectors),
np.asarray(eigenvectors) * np.asarray(eigenvalues)[..., None, :],
rtol=get_rtol(eigenvalues.dtype))
# 断言验证 matmul(a, eigenvectors) 几乎等于 eigenvectors * eigenvalues，考虑浮点数的相对误差
assert_(consistent_subclass(eigenvectors, a))
# 断言验证 eigenvectors 的类型与 a 一致
class TestEig(EigCases):
# 定义测试类 TestEig，继承自 EigCases
# 使用 pytest 的参数化装饰器，以便为每种数据类型运行此测试函数
@pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
def test_types(self, dtype):
# 创建一个二维数组 x，指定数据类型为 dtype
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
# 计算 x 的特征值和特征向量
w, v = np.linalg.eig(x)
# 断言特征值的数据类型与指定的 dtype 相同
assert_equal(w.dtype, dtype)
# 断言特征向量的数据类型与指定的 dtype 相同
assert_equal(v.dtype, dtype)
# 创建另一个二维数组 x，指定数据类型为 dtype
x = np.array([[1, 0.5], [-1, 1]], dtype=dtype)
# 计算 x 的特征值和特征向量
w, v = np.linalg.eig(x)
# 断言特征值的数据类型为由 get_complex_dtype 函数返回的复数数据类型
assert_equal(w.dtype, get_complex_dtype(dtype))
# 断言特征向量的数据类型为由 get_complex_dtype 函数返回的复数数据类型
assert_equal(v.dtype, get_complex_dtype(dtype))
def test_0_size(self):
# 检查各种零大小的数组是否正常工作
# 定义一个继承自 np.ndarray 的数组子类
class ArraySubclass(np.ndarray):
pass
# 创建一个形状为 (0, 1, 1)，数据类型为 np.int_ 的全零数组，并将其视图转换为 ArraySubclass 类型
a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
# 计算 a 的特征值和特征向量
res, res_v = linalg.eig(a)
# 断言特征向量的数据类型为 np.float64
assert_(res_v.dtype.type is np.float64)
# 断言特征值的数据类型为 np.float64
assert_(res.dtype.type is np.float64)
# 断言 a 的形状与 res_v 的形状相等
assert_equal(a.shape, res_v.shape)
# 断言 res 的形状为 (0, 1)
assert_equal((0, 1), res.shape)
# 这仅用于文档，可能需要更改：
# 断言 a 是 np.ndarray 类型的实例
assert_(isinstance(a, np.ndarray))
# 创建一个形状为 (0, 0)，数据类型为 np.complex64 的全零数组，并将其视图转换为 ArraySubclass 类型
a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
# 计算 a 的特征值和特征向量
res, res_v = linalg.eig(a)
# 断言特征向量的数据类型为 np.complex64
assert_(res_v.dtype.type is np.complex64)
# 断言特征值的数据类型为 np.complex64
assert_(res.dtype.type is np.complex64)
# 断言 a 的形状与 res_v 的形状相等
assert_equal(a.shape, res_v.shape)
# 断言 res 的形状为 (0,)
assert_equal((0,), res.shape)
# 这仅用于文档，可能需要更改：
# 断言 a 是 np.ndarray 类型的实例
assert_(isinstance(a, np.ndarray))
class SVDBaseTests:
hermitian = False
@pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
# 使用 pytest 的参数化装饰器，针对不同的 dtype 参数化测试类型
def test_types(self, dtype):
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
# 创建一个 numpy 数组 x，指定数据类型为 dtype
res = linalg.svd(x)
# 对 x 进行奇异值分解，并将结果存储在 res 中
U, S, Vh = res.U, res.S, res.Vh
# 从 res 中获取奇异值分解的结果 U、S、Vh
assert_equal(U.dtype, dtype)
# 断言 U 的数据类型与指定的 dtype 相同
assert_equal(S.dtype, get_real_dtype(dtype))
# 断言 S 的数据类型为 dtype 对应的实数类型
assert_equal(Vh.dtype, dtype)
# 断言 Vh 的数据类型与指定的 dtype 相同
s = linalg.svd(x, compute_uv=False, hermitian=self.hermitian)
# 对 x 进行奇异值分解，但不计算 U 和 Vh，同时考虑是否为共轭（当 hermitian=True 时）
assert_equal(s.dtype, get_real_dtype(dtype))
# 断言 s 的数据类型为 dtype 对应的实数类型
class SVDCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
def do(self, a, b, tags):
u, s, vt = linalg.svd(a, False)
# 对输入的矩阵 a 进行奇异值分解，不计算 U 和 V^T
assert_allclose(a, matmul(np.asarray(u) * np.asarray(s)[..., None, :],
np.asarray(vt)),
rtol=get_rtol(u.dtype))
# 断言 a 与通过 U、S、V^T 重构的矩阵相似
assert_(consistent_subclass(u, a))
# 断言 u 是 a 的一致子类
assert_(consistent_subclass(vt, a))
# 断言 vt 是 a 的一致子类
class TestSVD(SVDCases, SVDBaseTests):
def test_empty_identity(self):
""" Empty input should put an identity matrix in u or vh """
x = np.empty((4, 0))
# 创建一个空的 numpy 数组 x，形状为 (4, 0)
u, s, vh = linalg.svd(x, compute_uv=True, hermitian=self.hermitian)
# 对 x 进行奇异值分解，计算 U 和 V^H，同时考虑是否为共轭（当 hermitian=True 时）
assert_equal(u.shape, (4, 4))
# 断言 U 的形状为 (4, 4)
assert_equal(vh.shape, (0, 0))
# 断言 V^H 的形状为 (0, 0)
assert_equal(u, np.eye(4))
# 断言 U 与 4x4 单位矩阵相等
x = np.empty((0, 4))
# 创建一个空的 numpy 数组 x，形状为 (0, 4)
u, s, vh = linalg.svd(x, compute_uv=True, hermitian=self.hermitian)
# 对 x 进行奇异值分解，计算 U 和 V^H，同时考虑是否为共轭（当 hermitian=True 时）
assert_equal(u.shape, (0, 0))
# 断言 U 的形状为 (0, 0)
assert_equal(vh.shape, (4, 4))
# 断言 V^H 的形状为 (4, 4)
assert_equal(vh, np.eye(4))
# 断言 V^H 与 4x4 单位矩阵相等
def test_svdvals(self):
x = np.array([[1, 0.5], [0.5, 1]])
# 创建一个 numpy 数组 x
s_from_svd = linalg.svd(x, compute_uv=False, hermitian=self.hermitian)
# 对 x 进行奇异值分解，但不计算 U 和 V^H，同时考虑是否为共轭（当 hermitian=True 时）
s_from_svdvals = linalg.svdvals(x)
# 计算 x 的奇异值
assert_almost_equal(s_from_svd, s_from_svdvals)
# 断言通过奇异值分解和 svdvals 函数得到的奇异值近似相等
class SVDHermitianCases(HermitianTestCase, HermitianGeneralizedTestCase):
def do(self, a, b, tags):
u, s, vt = linalg.svd(a, False, hermitian=True)
# 对输入的矩阵 a 进行奇异值分解，不计算 U 和 V^H，并且假设为共轭
assert_allclose(a, matmul(np.asarray(u) * np.asarray(s)[..., None, :],
np.asarray(vt)),
rtol=get_rtol(u.dtype))
# 断言 a 与通过 U、S、V^H 重构的矩阵相似
def hermitian(mat):
axes = list(range(mat.ndim))
axes[-1], axes[-2] = axes[-2], axes[-1]
return np.conj(np.transpose(mat, axes=axes))
assert_almost_equal(np.matmul(u, hermitian(u)), np.broadcast_to(np.eye(u.shape[-1]), u.shape))
# 断言 U 与其共轭转置的乘积近似为单位矩阵
assert_almost_equal(np.matmul(vt, hermitian(vt)), np.broadcast_to(np.eye(vt.shape[-1]), vt.shape))
# 断言 V^H 与其共轭转置的乘积近似为单位矩阵
assert_equal(np.sort(s)[..., ::-1], s)
# 断言 s 为非增序列
assert_(consistent_subclass(u, a))
# 断言 u 是 a 的一致子类
assert_(consistent_subclass(vt, a))
# 断言 vt 是 a 的一致子类
class TestSVDHermitian(SVDHermitianCases, SVDBaseTests):
hermitian = True
class CondCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
# cond(x, p) for p in (None, 2, -2)
# 定义一个方法 `do`，接受三个参数 `a`, `b`, `tags`
def do(self, a, b, tags):
# 将参数 `a` 转换为 NumPy 数组 `c`，假设 `a` 可能是一个矩阵
c = asarray(a) # a might be a matrix
# 如果标签中包含 'size-0'
if 'size-0' in tags:
# 断言会抛出 LinAlgError 异常，测试矩阵 `c` 的条件数
assert_raises(LinAlgError, linalg.cond, c)
return
# 计算矩阵 `c` 的奇异值，不计算左右奇异向量
s = linalg.svd(c, compute_uv=False)
# 断言矩阵 `a` 的条件数近似等于最大和最小奇异值比值
assert_almost_equal(
linalg.cond(a), s[..., 0] / s[..., -1],
single_decimal=5, double_decimal=11)
# 断言矩阵 `a` 的 2-范数条件数近似等于最大和最小奇异值比值
assert_almost_equal(
linalg.cond(a, 2), s[..., 0] / s[..., -1],
single_decimal=5, double_decimal=11)
# 断言矩阵 `a` 的 -2-范数条件数近似等于最小和最大奇异值比值的倒数
assert_almost_equal(
linalg.cond(a, -2), s[..., -1] / s[..., 0],
single_decimal=5, double_decimal=11)
# 计算矩阵 `c` 的逆矩阵
cinv = np.linalg.inv(c)
# 断言矩阵 `a` 的 1-范数条件数近似等于绝对值矩阵 `c` 按行求和后的最大值乘以逆矩阵的相同计算结果
assert_almost_equal(
linalg.cond(a, 1),
abs(c).sum(-2).max(-1) * abs(cinv).sum(-2).max(-1),
single_decimal=5, double_decimal=11)
# 断言矩阵 `a` 的 -1-范数条件数近似等于绝对值矩阵 `c` 按行求和后的最小值乘以逆矩阵的相同计算结果
assert_almost_equal(
linalg.cond(a, -1),
abs(c).sum(-2).min(-1) * abs(cinv).sum(-2).min(-1),
single_decimal=5, double_decimal=11)
# 断言矩阵 `a` 的无穷范数条件数近似等于绝对值矩阵 `c` 按列求和后的最大值乘以逆矩阵的相同计算结果
assert_almost_equal(
linalg.cond(a, np.inf),
abs(c).sum(-1).max(-1) * abs(cinv).sum(-1).max(-1),
single_decimal=5, double_decimal=11)
# 断言矩阵 `a` 的负无穷范数条件数近似等于绝对值矩阵 `c` 按列求和后的最小值乘以逆矩阵的相同计算结果
assert_almost_equal(
linalg.cond(a, -np.inf),
abs(c).sum(-1).min(-1) * abs(cinv).sum(-1).min(-1),
single_decimal=5, double_decimal=11)
# 断言矩阵 `a` 的 Frobenius 范数条件数近似等于矩阵 `c` 的 Frobenius 范数平方和乘以逆矩阵的相同计算结果的平方根
assert_almost_equal(
linalg.cond(a, 'fro'),
np.sqrt((abs(c)**2).sum(-1).sum(-1)
* (abs(cinv)**2).sum(-1).sum(-1)),
single_decimal=5, double_decimal=11)
# 定义一个名为 TestCond 的类，继承自 CondCases
class TestCond(CondCases):
# 定义一个名为 test_basic_nonsvd 的方法
def test_basic_nonsvd(self):
# 对非奇异值分解进行烟雾测试
A = array([[1., 0, 1], [0, -2., 0], [0, 0, 3.]])
# 使用无穷范数计算矩阵 A 的条件数，并与 4 进行差值比较
assert_almost_equal(linalg.cond(A, inf), 4)
# 使用负无穷范数计算矩阵 A 的条件数，并与 2/3 进行差值比较
assert_almost_equal(linalg.cond(A, -inf), 2/3)
# 使用 1 范数计算矩阵 A 的条件数，并与 4 进行差值比较
assert_almost_equal(linalg.cond(A, 1), 4)
# 使用 -1 范数计算矩阵 A 的条件数，并与 0.5 进行差值比较
assert_almost_equal(linalg.cond(A, -1), 0.5)
# 使用 Frobenius 范数计算矩阵 A 的条件数，并与 sqrt(265 / 12) 进行差值比较
assert_almost_equal(linalg.cond(A, 'fro'), np.sqrt(265 / 12))
# 定义一个名为 test_singular 的方法
def test_singular(self):
# 奇异矩阵对于正范数有无限条件数，而负范数不应触发异常
As = [np.zeros((2, 2)), np.ones((2, 2))]
p_pos = [None, 1, 2, 'fro']
p_neg = [-1, -2]
for A, p in itertools.product(As, p_pos):
# 反转可能不会达到确切的无穷大，因此只需检查数值非常大
assert_(linalg.cond(A, p) > 1e15)
for A, p in itertools.product(As, p_neg):
linalg.cond(A, p)
# 定义一个名为 test_nan 的方法
@pytest.mark.xfail(True, run=False,
reason="Platform/LAPACK-dependent failure, "
"see gh-18914")
def test_nan(self):
# nan 应该被传递，而不是转换为无穷大
ps = [None, 1, -1, 2, -2, 'fro']
p_pos = [None, 1, 2, 'fro']
A = np.ones((2, 2))
A[0,1] = np.nan
for p in ps:
c = linalg.cond(A, p)
assert_(isinstance(c, np.float64))
assert_(np.isnan(c))
A = np.ones((3, 2, 2))
A[1,0,1] = np.nan
for p in ps:
c = linalg.cond(A, p)
assert_(np.isnan(c[1]))
if p in p_pos:
assert_(c[0] > 1e15)
assert_(c[2] > 1e15)
else:
assert_(not np.isnan(c[0]))
assert_(not np.isnan(c[2]))
# 定义一个名为 test_stacked_singular 的方法
def test_stacked_singular(self):
# 当只有部分堆叠矩阵是奇异时，检查行为
np.random.seed(1234)
A = np.random.rand(2, 2, 2, 2)
A[0,0] = 0
A[1,1] = 0
for p in (None, 1, 2, 'fro', -1, -2):
c = linalg.cond(A, p)
assert_equal(c[0,0], np.inf)
assert_equal(c[1,1], np.inf)
assert_(np.isfinite(c[0,1]))
assert_(np.isfinite(c[1,0]))
# 定义一个名为 PinvCases 的类，继承自多个测试用例类
class PinvCases(LinalgSquareTestCase,
LinalgNonsquareTestCase,
LinalgGeneralizedSquareTestCase,
LinalgGeneralizedNonsquareTestCase):
# 定义一个 do 方法，接受 a, b, tags 作为参数
def do(self, a, b, tags):
# 计算矩阵 a 的伪逆
a_ginv = linalg.pinv(a)
# `a @ a_ginv == I` 如果 a 是奇异的话可能不成立
dot = matmul
# 断言矩阵乘积 a @ a_ginv 与 a 之间的关系
assert_almost_equal(dot(dot(a, a_ginv), a), a, single_decimal=5, double_decimal=11)
# 断言 a_ginv 与 a 具有一致的子类
assert_(consistent_subclass(a_ginv, a))
# 定义一个名为 TestPinv 的类，继承自 PinvCases
class TestPinv(PinvCases):
pass
# 定义一个名为 PinvHermitianCases 的类，继承自 HermitianTestCase 和 HermitianGeneralizedTestCase
class PinvHermitianCases(HermitianTestCase, HermitianGeneralizedTestCase):
def do(self, a, b, tags):
# 计算矩阵 a 的广义逆
a_ginv = linalg.pinv(a, hermitian=True)
# 检查 `a @ a_ginv == I` 是否成立，如果矩阵 a 是奇异的，则不成立
dot = matmul
# 断言验证 dot(dot(a, a_ginv), a) 接近于 a，精确到小数点后五位和十一位
assert_almost_equal(dot(dot(a, a_ginv), a), a, single_decimal=5, double_decimal=11)
# 断言验证 a_ginv 是否是 a 的一致子类
assert_(consistent_subclass(a_ginv, a))
class TestPinvHermitian(PinvHermitianCases):
pass
def test_pinv_rtol_arg():
# 创建一个测试用的二维数组 `a`
a = np.array([[1, 2, 3], [4, 1, 1], [2, 3, 1]])
# 使用 `rcond` 参数计算广义逆矩阵，并使用 `rtol` 参数计算广义逆矩阵，比较它们的近似值
assert_almost_equal(
np.linalg.pinv(a, rcond=0.5),
np.linalg.pinv(a, rtol=0.5),
)
# 当 `rtol` 和 `rcond` 同时设置时，使用 `pytest.raises` 检查是否会引发 ValueError 异常
with pytest.raises(
ValueError, match=r"`rtol` and `rcond` can't be both set."
):
np.linalg.pinv(a, rcond=0.5, rtol=0.5)
class DetCases(LinalgSquareTestCase, LinalgGeneralizedSquareTestCase):
def do(self, a, b, tags):
# 计算矩阵的行列式值 `d`
d = linalg.det(a)
# 计算矩阵的行列式的符号和对数绝对值
res = linalg.slogdet(a)
s, ld = res.sign, res.logabsdet
# 根据矩阵的数据类型，转换为适当的复数类型
if asarray(a).dtype.type in (single, double):
ad = asarray(a).astype(double)
else:
ad = asarray(a).astype(cdouble)
# 计算矩阵的特征值
ev = linalg.eigvals(ad)
# 检查行列式的计算结果与特征值乘积的一致性
assert_almost_equal(d, multiply.reduce(ev, axis=-1))
# 检查行列式符号乘以指数对数绝对值是否与特征值乘积的一致性
assert_almost_equal(s * np.exp(ld), multiply.reduce(ev, axis=-1))
# 将符号 `s` 和对数绝对值 `ld` 至少转换为一维数组
s = np.atleast_1d(s)
ld = np.atleast_1d(ld)
# 根据符号 `s` 的非零判断条件，检查其绝对值是否为 1
m = (s != 0)
assert_almost_equal(np.abs(s[m]), 1)
# 对于非零符号 `s`，检查对数绝对值 `ld` 是否为负无穷
assert_equal(ld[~m], -inf)
class TestDet(DetCases):
def test_zero(self):
# 对于零矩阵的行列式和特殊情况进行检查
assert_equal(linalg.det([[0.0]]), 0.0)
assert_equal(type(linalg.det([[0.0]])), double)
assert_equal(linalg.det([[0.0j]]), 0.0)
assert_equal(type(linalg.det([[0.0j]])), cdouble)
assert_equal(linalg.slogdet([[0.0]]), (0.0, -inf))
assert_equal(type(linalg.slogdet([[0.0]])[0]), double)
assert_equal(type(linalg.slogdet([[0.0]])[1]), double)
assert_equal(linalg.slogdet([[0.0j]]), (0.0j, -inf))
assert_equal(type(linalg.slogdet([[0.0j]])[0]), cdouble)
assert_equal(type(linalg.slogdet([[0.0j]])[1]), double)
@pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
def test_types(self, dtype):
# 对于不同数据类型的矩阵进行行列式和特征值计算的数据类型检查
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
assert_equal(np.linalg.det(x).dtype, dtype)
ph, s = np.linalg.slogdet(x)
assert_equal(s.dtype, get_real_dtype(dtype))
assert_equal(ph.dtype, dtype)
def test_0_size(self):
# 对于大小为 0 的矩阵，行列式和特征值计算的特殊情况检查
a = np.zeros((0, 0), dtype=np.complex64)
res = linalg.det(a)
assert_equal(res, 1.)
assert_(res.dtype.type is np.complex64)
res = linalg.slogdet(a)
assert_equal(res, (1, 0))
assert_(res[0].dtype.type is np.complex64)
assert_(res[1].dtype.type is np.float32)
a = np.zeros((0, 0), dtype=np.float64)
res = linalg.det(a)
assert_equal(res, 1.)
assert_(res.dtype.type is np.float64)
res = linalg.slogdet(a)
assert_equal(res, (1, 0))
assert_(res[0].dtype.type is np.float64)
assert_(res[1].dtype.type is np.float64)
# 定义一个方法，接受三个参数：a，b，tags
def do(self, a, b, tags):
# 将a转换为NumPy数组
arr = np.asarray(a)
# 获取数组的形状 m, n
m, n = arr.shape
# 对a进行奇异值分解，返回结果是u, s, vt
u, s, vt = linalg.svd(a, False)
# 使用最小二乘法求解线性方程组a*x=b，返回结果是x, residuals, rank, sv
x, residuals, rank, sv = linalg.lstsq(a, b, rcond=-1)
# 如果m等于0
if m == 0:
# 断言x的所有元素都为0
assert_((x == 0).all())
# 如果m小于等于n
if m <= n:
# 断言b与a*x的点积近似相等
assert_almost_equal(b, dot(a, x))
# 断言rank等于m
assert_equal(rank, m)
else:
# 断言rank等于n
assert_equal(rank, n)
# 断言sv近似等于s
assert_almost_equal(sv, sv.__array_wrap__(s))
# 如果rank等于n且m大于n
if rank == n and m > n:
# 预期残差的平方和
expect_resids = (
np.asarray(abs(np.dot(a, x) - b)) ** 2).sum(axis=0)
# 将expect_resids转换为NumPy数组
expect_resids = np.asarray(expect_resids)
# 如果b是一维数组
if np.asarray(b).ndim == 1:
# 调整expect_resids的形状
expect_resids.shape = (1,)
# 断言residuals与expect_resids的形状相同
assert_equal(residuals.shape, expect_resids.shape)
else:
# 创建一个空的NumPy数组，其数据类型与x相同
expect_resids = np.array([]).view(type(x))
# 断言residuals近似等于expect_resids
assert_almost_equal(residuals, expect_resids)
# 断言residuals的数据类型为浮点数
assert_(np.issubdtype(residuals.dtype, np.floating))
# 断言x和b具有一致的子类类型
assert_(consistent_subclass(x, b))
# 断言residuals和b具有一致的子类类型
assert_(consistent_subclass(residuals, b))
class TestLstsq(LstsqCases):
# 继承自 LstsqCases 类的测试类 TestLstsq
def test_rcond(self):
# 测试函数，用于测试 linalg.lstsq 函数的 rcond 参数
a = np.array([[0., 1., 0., 1., 2., 0.],
[0., 2., 0., 0., 1., 0.],
[1., 0., 1., 0., 0., 4.],
[0., 0., 0., 2., 3., 0.]]).T
# 创建一个 numpy 数组 a，表示系数矩阵的转置
b = np.array([1, 0, 0, 0, 0, 0])
# 创建一个 numpy 数组 b，表示常数向量
x, residuals, rank, s = linalg.lstsq(a, b, rcond=-1)
# 使用 linalg.lstsq 函数求解线性方程组，rcond 参数为 -1
assert_(rank == 4)
# 断言 rank 等于 4
x, residuals, rank, s = linalg.lstsq(a, b)
# 再次使用 linalg.lstsq 函数求解线性方程组，默认 rcond 参数
assert_(rank == 3)
# 断言 rank 等于 3
x, residuals, rank, s = linalg.lstsq(a, b, rcond=None)
# 第三次使用 linalg.lstsq 函数求解线性方程组，rcond 参数为 None
assert_(rank == 3)
# 断言 rank 等于 3
@pytest.mark.parametrize(["m", "n", "n_rhs"], [
(4, 2, 2),
(0, 4, 1),
(0, 4, 2),
(4, 0, 1),
(4, 0, 2),
(4, 2, 0),
(0, 0, 0)
])
def test_empty_a_b(self, m, n, n_rhs):
# 参数化测试函数，用于测试当 a 或 b 为空时的情况
a = np.arange(m * n).reshape(m, n)
# 创建一个 m 行 n 列的 numpy 数组 a
b = np.ones((m, n_rhs))
# 创建一个元素全为 1 的 m 行 n_rhs 列的 numpy 数组 b
x, residuals, rank, s = linalg.lstsq(a, b, rcond=None)
# 使用 linalg.lstsq 函数求解线性方程组，rcond 参数为 None
if m == 0:
assert_((x == 0).all())
# 如果 m 为 0，则断言 x 中所有元素为 0
assert_equal(x.shape, (n, n_rhs))
# 断言 x 的形状为 (n, n_rhs)
assert_equal(residuals.shape, ((n_rhs,) if m > n else (0,)))
# 断言 residuals 的形状符合预期 ((n_rhs,) 或 (0,))
if m > n and n_rhs > 0:
# 如果 m 大于 n 并且 n_rhs 大于 0
# residuals 正好是 b 列的平方范数
r = b - np.dot(a, x)
assert_almost_equal(residuals, (r * r).sum(axis=-2))
# 断言 residuals 与 (r * r).sum(axis=-2) 几乎相等
assert_equal(rank, min(m, n))
# 断言 rank 等于 m 和 n 中的较小值
assert_equal(s.shape, (min(m, n),))
# 断言 s 的形状为 (min(m, n),)
def test_incompatible_dims(self):
# 测试函数，用于测试当维度不兼容时的情况
x = np.array([0, 1, 2, 3])
# 创建一个 numpy 数组 x
y = np.array([-1, 0.2, 0.9, 2.1, 3.3])
# 创建一个 numpy 数组 y
A = np.vstack([x, np.ones(len(x))]).T
# 创建一个由 x 和一个全为 1 的数组堆叠而成的 numpy 数组 A
with assert_raises_regex(LinAlgError, "Incompatible dimensions"):
linalg.lstsq(A, y, rcond=None)
# 使用 linalg.lstsq 函数求解线性方程组，rcond 参数为 None
@pytest.mark.parametrize('dt', [np.dtype(c) for c in '?bBhHiIqQefdgFDGO'])
class TestMatrixPower:
# 参数化测试类，用于测试 matrix_power 函数的不同数据类型
rshft_0 = np.eye(4)
# 创建一个 4x4 的单位矩阵赋值给 rshft_0
rshft_1 = rshft_0[[3, 0, 1, 2]]
# 创建 rshft_0 的第 4 行换到第 1 行的矩阵赋值给 rshft_1
rshft_2 = rshft_0[[2, 3, 0, 1]]
# 创建 rshft_0 的第 3 行换到第 0 行的矩阵赋值给 rshft_2
rshft_3 = rshft_0[[1, 2, 3, 0]]
# 创建 rshft_0 的第 2 行换到第 3 行的矩阵赋值给 rshft_3
rshft_all = [rshft_0, rshft_1, rshft_2, rshft_3]
# 创建包含 rshft_0、rshft_1、rshft_2 和 rshft_3 的列表赋值给 rshft_all
noninv = array([[1, 0], [0, 0]])
# 创建一个非逆矩阵赋值给 noninv
stacked = np.block([[[rshft_0]]]*2)
# 创建一个堆叠 rshft_0 两次的 numpy 数组赋值给 stacked
# FIXME the 'e' dtype might work in future
dtnoinv = [object, np.dtype('e'), np.dtype('g'), np.dtype('G')]
# 创建包含 object、'e'、'g' 和 'G' 数据类型的列表赋值给 dtnoinv
def test_large_power(self, dt):
# 测试函数，用于测试 matrix_power 函数对于大数幂次的计算
rshft = self.rshft_1.astype(dt)
# 将 rshft_1 转换为指定数据类型 dt
assert_equal(
matrix_power(rshft, 2**100 + 2**10 + 2**5 + 0), self.rshft_0)
# 断言 matrix_power 计算结果与 rshft_0 在指定幂次下的相等性
assert_equal(
matrix_power(rshft, 2**100 + 2**10 + 2**5 + 1), self.rshft_1)
# 断言 matrix_power 计算结果与 rshft_1 在指定幂次下的相等性
assert_equal(
matrix_power(rshft, 2**100 + 2**10 + 2**5 + 2), self.rshft_2)
# 断言 matrix_power 计算结果与 rshft_2 在指定幂次下的相等性
assert_equal(
matrix_power(rshft, 2**100 + 2**10 + 2**5 + 3), self.rshft_3)
# 断言 matrix_power 计算结果与 rshft_3 在指定幂次下的相等性
def test_power_is_zero(self, dt):
# 测试函数，用于测试 matrix_power 函数对于幂次为 0 的情况
def tz(M):
mz = matrix_power(M, 0)
# 计算矩阵 M 的零次幂
assert_equal(mz, identity_like_generalized(M))
# 断言 mz 等于 M 的广义单位矩阵
assert_equal(mz.dtype, M.dtype)
# 断言 mz
# 定义一个测试方法，用于测试矩阵乘幂为1的情况
def test_power_is_one(self, dt):
# 定义内部函数 tz，接受一个矩阵 mat 作为参数
def tz(mat):
# 计算矩阵 mat 的一次幂
mz = matrix_power(mat, 1)
# 断言计算结果 mz 等于原始矩阵 mat
assert_equal(mz, mat)
# 断言计算结果 mz 的数据类型等于原始矩阵 mat 的数据类型
assert_equal(mz.dtype, mat.dtype)
# 遍历 self.rshft_all 中的每一个矩阵
for mat in self.rshft_all:
# 调用 tz 函数，传入 mat 转换为指定数据类型 dt 后的结果
tz(mat.astype(dt))
# 如果数据类型 dt 不是对象类型
if dt != object:
# 调用 tz 函数，传入 self.stacked 转换为指定数据类型 dt 后的结果
tz(self.stacked.astype(dt))
# 定义一个测试方法，用于测试矩阵乘幂为2的情况
def test_power_is_two(self, dt):
# 定义内部函数 tz，接受一个矩阵 mat 作为参数
def tz(mat):
# 计算矩阵 mat 的二次幂
mz = matrix_power(mat, 2)
# 根据矩阵 mat 的数据类型选择使用 matmul 或 dot 函数进行矩阵乘法
mmul = matmul if mat.dtype != object else dot
# 断言计算结果 mz 等于 mat 与自身的乘积
assert_equal(mz, mmul(mat, mat))
# 断言计算结果 mz 的数据类型等于原始矩阵 mat 的数据类型
assert_equal(mz.dtype, mat.dtype)
# 遍历 self.rshft_all 中的每一个矩阵
for mat in self.rshft_all:
# 调用 tz 函数，传入 mat 转换为指定数据类型 dt 后的结果
tz(mat.astype(dt))
# 如果数据类型 dt 不是对象类型
if dt != object:
# 调用 tz 函数，传入 self.stacked 转换为指定数据类型 dt 后的结果
tz(self.stacked.astype(dt))
# 定义一个测试方法，用于测试矩阵乘幂为-1的情况
def test_power_is_minus_one(self, dt):
# 定义内部函数 tz，接受一个矩阵 mat 作为参数
def tz(mat):
# 计算矩阵 mat 的逆矩阵
invmat = matrix_power(mat, -1)
# 根据矩阵 mat 的数据类型选择使用 matmul 或 dot 函数进行矩阵乘法
mmul = matmul if mat.dtype != object else dot
# 断言计算结果 invmat 与 mat 的乘积接近单位矩阵
assert_almost_equal(
mmul(invmat, mat), identity_like_generalized(mat))
# 遍历 self.rshft_all 中的每一个矩阵
for mat in self.rshft_all:
# 如果数据类型 dt 不在 self.dtnoinv 中
if dt not in self.dtnoinv:
# 调用 tz 函数，传入 mat 转换为指定数据类型 dt 后的结果
tz(mat.astype(dt))
# 定义一个测试方法，用于测试非法的乘幂情况
def test_exceptions_bad_power(self, dt):
# 将 self.rshft_0 转换为指定数据类型 dt 后赋值给 mat
mat = self.rshft_0.astype(dt)
# 断言对于非整数乘幂会抛出 TypeError 异常
assert_raises(TypeError, matrix_power, mat, 1.5)
assert_raises(TypeError, matrix_power, mat, [1])
# 定义一个测试方法，用于测试非方阵矩阵的乘幂情况
def test_exceptions_non_square(self, dt):
# 断言对于非方阵矩阵会抛出 LinAlgError 异常
assert_raises(LinAlgError, matrix_power, np.array([1], dt), 1)
assert_raises(LinAlgError, matrix_power, np.array([[1], [2]], dt), 1)
assert_raises(LinAlgError, matrix_power, np.ones((4, 3, 2), dt), 1)
# 使用 pytest.mark.skipif 装饰器标记一个测试方法，用于测试不可逆矩阵的乘幂情况
@pytest.mark.skipif(IS_WASM, reason="fp errors don't work in wasm")
def test_exceptions_not_invertible(self, dt):
# 如果数据类型 dt 在 self.dtnoinv 中，则直接返回
if dt in self.dtnoinv:
return
# 将 self.noninv 转换为指定数据类型 dt 后赋值给 mat
mat = self.noninv.astype(dt)
# 断言对于不可逆矩阵会抛出 LinAlgError 异常
assert_raises(LinAlgError, matrix_power, mat, -1)
class TestEigvalshCases(HermitianTestCase, HermitianGeneralizedTestCase):
def do(self, a, b, tags):
# note that eigenvalue arrays returned by eig must be sorted since
# their order isn't guaranteed.
# 使用 np.linalg.eigvalsh 计算矩阵 a 的特征值，并要求以升序排列
ev = linalg.eigvalsh(a, 'L')
# 使用 np.linalg.eig 计算矩阵 a 的特征值和特征向量
evalues, evectors = linalg.eig(a)
# 对通过 np.linalg.eig 计算得到的特征值进行排序
evalues.sort(axis=-1)
# 断言 np.linalg.eigvalsh 计算得到的特征值数组与 np.linalg.eig 计算得到的特征值数组近似相等
assert_allclose(ev, evalues, rtol=get_rtol(ev.dtype))
# 使用 np.linalg.eigvalsh 计算矩阵 a 的特征值，并要求以降序排列
ev2 = linalg.eigvalsh(a, 'U')
# 断言 np.linalg.eigvalsh 计算得到的特征值数组与先前计算的 evalues 近似相等
assert_allclose(ev2, evalues, rtol=get_rtol(ev.dtype))
class TestEigvalsh:
@pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
def test_types(self, dtype):
# 创建一个指定数据类型的数组 x
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
# 使用 np.linalg.eigvalsh 计算数组 x 的特征值
w = np.linalg.eigvalsh(x)
# 断言计算得到的特征值数组的数据类型与预期的实数类型一致
assert_equal(w.dtype, get_real_dtype(dtype))
def test_invalid(self):
# 创建一个指定数据类型的数组 x
x = np.array([[1, 0.5], [0.5, 1]], dtype=np.float32)
# 断言当使用无效的 UPLO 参数 'lrong' 调用 np.linalg.eigvalsh 时会引发 ValueError
assert_raises(ValueError, np.linalg.eigvalsh, x, UPLO="lrong")
# 断言当使用 UPLO 参数为 "lower" 调用 np.linalg.eigvalsh 时会引发 ValueError
assert_raises(ValueError, np.linalg.eigvalsh, x, "lower")
# 断言当使用 UPLO 参数为 "upper" 调用 np.linalg.eigvalsh 时会引发 ValueError
assert_raises(ValueError, np.linalg.eigvalsh, x, "upper")
def test_UPLO(self):
# 创建两个指定数据类型的矩阵 Klo 和 Kup，分别用于测试 'L' 和 'U' 选项
Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
# 创建一个目标特征值数组 tgt 和比较的相对误差 rtol
tgt = np.array([-1, 1], dtype=np.double)
rtol = get_rtol(np.double)
# 检查默认情况下 UPLO 参数为 'L' 时 np.linalg.eigvalsh 返回的特征值与目标 tgt 近似相等
w = np.linalg.eigvalsh(Klo)
assert_allclose(w, tgt, rtol=rtol)
# 检查 UPLO 参数为 'L' 时 np.linalg.eigvalsh 返回的特征值与目标 tgt 近似相等
w = np.linalg.eigvalsh(Klo, UPLO='L')
assert_allclose(w, tgt, rtol=rtol)
# 检查 UPLO 参数为 'l' 时 np.linalg.eigvalsh 返回的特征值与目标 tgt 近似相等
w = np.linalg.eigvalsh(Klo, UPLO='l')
assert_allclose(w, tgt, rtol=rtol)
# 检查 UPLO 参数为 'U' 时 np.linalg.eigvalsh 返回的特征值与目标 tgt 近似相等
w = np.linalg.eigvalsh(Kup, UPLO='U')
assert_allclose(w, tgt, rtol=rtol)
# 检查 UPLO 参数为 'u' 时 np.linalg.eigvalsh 返回的特征值与目标 tgt 近似相等
w = np.linalg.eigvalsh(Kup, UPLO='u')
assert_allclose(w, tgt, rtol=rtol)
def test_0_size(self):
# 检查所有大小为 0 的数组类型都能正常工作
# 创建一个 ArraySubclass 类型的零大小数组 a，并将其转换为 np.ndarray 类型
class ArraySubclass(np.ndarray):
pass
a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
# 使用 np.linalg.eigvalsh 计算零大小数组 a 的特征值
res = linalg.eigvalsh(a)
# 断言计算得到的特征值数组的数据类型为 np.float64
assert_(res.dtype.type is np.float64)
# 断言计算得到的特征值数组的形状为 (0, 1)
assert_equal((0, 1), res.shape)
# 这只是为了说明，可能需要更改：
# 断言计算得到的 res 是 np.ndarray 类型的实例
assert_(isinstance(res, np.ndarray))
a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
# 使用 np.linalg.eigvalsh 计算零大小数组 a 的特征值
res = linalg.eigvalsh(a)
# 断言计算得到的特征值数组的数据类型为 np.float32
assert_(res.dtype.type is np.float32)
# 断言计算得到的特征值数组的形状为 (0,)
assert_equal((0,), res.shape)
# 这只是为了说明，可能需要更改：
# 断言计算得到的 res 是 np.ndarray 类型的实例
assert_(isinstance(res, np.ndarray))
class TestEighCases(HermitianTestCase, HermitianGeneralizedTestCase):
# 未完待续，注释继续在下一个类的注释中添加
# 定义一个方法 `do`，接受参数 `a`, `b`, `tags`
def do(self, a, b, tags):
# 警告：eig 返回的特征值数组需要排序，因为它们的顺序不保证一致。
# 使用 linalg.eigh 计算矩阵 a 的特征值和特征向量
res = linalg.eigh(a)
# 分别获取特征值和特征向量
ev, evc = res.eigenvalues, res.eigenvectors
# 使用 linalg.eig 计算矩阵 a 的特征值和特征向量
evalues, evectors = linalg.eig(a)
# 对特征值进行排序，沿着最后一个轴进行排序
evalues.sort(axis=-1)
# 断言 ev（通过 linalg.eigh 得到的特征值）与 evalues（通过 linalg.eig 得到并排序的特征值）几乎相等
assert_almost_equal(ev, evalues)
# 断言矩阵 a 与特征向量 evc 的乘积，等于 evc 的每个列向量乘以对应的特征值，这里使用 matmul 进行矩阵乘法
assert_allclose(matmul(a, evc),
np.asarray(ev)[..., None, :] * np.asarray(evc),
rtol=get_rtol(ev.dtype))
# 使用 linalg.eigh 计算矩阵 a 的特征值和特征向量，参数 'U' 表示返回的特征向量是未正交化的
ev2, evc2 = linalg.eigh(a, 'U')
# 断言 ev2（通过 linalg.eigh 得到的未正交化的特征值）与 evalues（通过 linalg.eig 得到并排序的特征值）几乎相等
assert_almost_equal(ev2, evalues)
# 断言矩阵 a 与特征向量 evc2 的乘积，等于 evc2 的每个列向量乘以对应的特征值，这里使用 matmul 进行矩阵乘法
assert_allclose(matmul(a, evc2),
np.asarray(ev2)[..., None, :] * np.asarray(evc2),
rtol=get_rtol(ev.dtype), err_msg=repr(a))
class TestEigh:
# 使用 pytest 的参数化装饰器，测试不同的数据类型
@pytest.mark.parametrize('dtype', [single, double, csingle, cdouble])
def test_types(self, dtype):
# 创建一个二维数组 x，指定数据类型为 dtype
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
# 计算对称矩阵 x 的特征值和特征向量
w, v = np.linalg.eigh(x)
# 断言特征值的数据类型与给定 dtype 对应的实数部分的数据类型相同
assert_equal(w.dtype, get_real_dtype(dtype))
# 断言特征向量的数据类型与给定 dtype 相同
assert_equal(v.dtype, dtype)
# 测试处理无效输入的情况
def test_invalid(self):
# 创建一个二维数组 x，指定数据类型为 np.float32
x = np.array([[1, 0.5], [0.5, 1]], dtype=np.float32)
# 断言当 UPLO 参数为 "lrong" 时，调用 np.linalg.eigh(x) 抛出 ValueError 异常
assert_raises(ValueError, np.linalg.eigh, x, UPLO="lrong")
# 断言当 UPLO 参数为 "lower" 时，调用 np.linalg.eigh(x) 抛出 ValueError 异常
assert_raises(ValueError, np.linalg.eigh, x, "lower")
# 断言当 UPLO 参数为 "upper" 时，调用 np.linalg.eigh(x) 抛出 ValueError 异常
assert_raises(ValueError, np.linalg.eigh, x, "upper")
# 测试不同的 UPLO 参数对结果的影响
def test_UPLO(self):
# 创建两个双精度二维数组 Klo 和 Kup
Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
# 创建目标特征值数组 tgt，数据类型为双精度
tgt = np.array([-1, 1], dtype=np.double)
# 获取双精度数值的相对误差容限
rtol = get_rtol(np.double)
# 检查默认情况下 UPLO 参数为 'L' 时的特征值计算结果
w, v = np.linalg.eigh(Klo)
assert_allclose(w, tgt, rtol=rtol)
# 检查显式指定 UPLO 参数为 'L' 时的特征值计算结果
w, v = np.linalg.eigh(Klo, UPLO='L')
assert_allclose(w, tgt, rtol=rtol)
# 检查 UPLO 参数为 'l' 时的特征值计算结果
w, v = np.linalg.eigh(Klo, UPLO='l')
assert_allclose(w, tgt, rtol=rtol)
# 检查 UPLO 参数为 'U' 时的特征值计算结果
w, v = np.linalg.eigh(Kup, UPLO='U')
assert_allclose(w, tgt, rtol=rtol)
# 检查 UPLO 参数为 'u' 时的特征值计算结果
w, v = np.linalg.eigh(Kup, UPLO='u')
assert_allclose(w, tgt, rtol=rtol)
# 测试处理大小为零的情况
def test_0_size(self):
# 检查各种零大小数组的特征值和特征向量计算
class ArraySubclass(np.ndarray):
pass
# 创建一个 shape 为 (0, 1, 1) 的整数类型的零数组 a，视图类型为 ArraySubclass
a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
res, res_v = linalg.eigh(a)
# 断言特征向量的数据类型为 np.float64
assert_(res_v.dtype.type is np.float64)
# 断言特征值的数据类型为 np.float64
assert_(res.dtype.type is np.float64)
# 断言数组 a 的形状与特征向量 res_v 的形状相同
assert_equal(a.shape, res_v.shape)
# 断言特征值 res 的形状为 (0, 1)
assert_equal((0, 1), res.shape)
# 这是为了文档说明，可能会有变化：
assert_(isinstance(a, np.ndarray))
# 创建一个 shape 为 (0, 0) 的复数类型的零数组 a，视图类型为 ArraySubclass
a = np.zeros((0, 0), dtype=np.complex64).view(ArraySubclass)
res, res_v = linalg.eigh(a)
# 断言特征向量的数据类型为 np.complex64
assert_(res_v.dtype.type is np.complex64)
# 断言特征值的数据类型为 np.float32
assert_(res.dtype.type is np.float32)
# 断言数组 a 的形状与特征向量 res_v 的形状相同
assert_equal(a.shape, res_v.shape)
# 断言特征值 res 的形状为 (0,)
assert_equal((0,), res.shape)
# 这是为了文档说明，可能会有变化：
assert_(isinstance(a, np.ndarray))
class _TestNormBase:
dt = None
dec = None
# 静态方法：检查输入数组的数据类型
@staticmethod
def check_dtype(x, res):
# 如果输入数组 x 的数据类型是 np.inexact 的子类
if issubclass(x.dtype.type, np.inexact):
# 断言结果数组 res 的数据类型与 x 实数部分的数据类型相同
assert_equal(res.dtype, x.real.dtype)
else:
# 对于整数输入，不必测试输出的浮点精度。
assert_(issubclass(res.dtype.type, np.floating))
class _TestNormGeneral(_TestNormBase):
# 测试空数组的情况
def test_empty(self):
# 断言对空列表计算的范数为 0.0
assert_equal(norm([]), 0.0)
# 断言对 shape 为空的数组计算的范数为 0.0
assert_equal(norm(array([], dtype=self.dt)), 0.0)
# 断言对至少二维 shape 为空的数组计算的范数为 0.0
assert_equal(norm(atleast_2d(array([], dtype=self.dt))), 0.0)
# 定义测试向量返回类型的方法
def test_vector_return_type(self):
# 创建一个 NumPy 数组 [1, 0, 1]
a = np.array([1, 0, 1])
# 获取所有精确类型的类型码
exact_types = np.typecodes['AllInteger']
# 获取所有非精确类型的类型码
inexact_types = np.typecodes['AllFloat']
# 合并精确和非精确类型的类型码
all_types = exact_types + inexact_types
# 遍历所有类型
for each_type in all_types:
# 将数组 a 转换为当前遍历的类型
at = a.astype(each_type)
# 计算 L-infinity 范数
an = norm(at, -np.inf)
self.check_dtype(at, an)
# 断言 L-infinity 范数近似为 0.0
assert_almost_equal(an, 0.0)
# 忽略运行时警告并计算 L-1 范数
with suppress_warnings() as sup:
sup.filter(RuntimeWarning, "divide by zero encountered")
an = norm(at, -1)
self.check_dtype(at, an)
# 断言 L-1 范数近似为 0.0
assert_almost_equal(an, 0.0)
# 计算 L-0 范数
an = norm(at, 0)
self.check_dtype(at, an)
# 断言 L-0 范数为 2
assert_almost_equal(an, 2)
# 计算 L-1 范数
an = norm(at, 1)
self.check_dtype(at, an)
# 断言 L-1 范数为 2.0
assert_almost_equal(an, 2.0)
# 计算 L-2 范数
an = norm(at, 2)
self.check_dtype(at, an)
# 断言 L-2 范数近似为 sqrt(2)
assert_almost_equal(an, an.dtype.type(2.0)**an.dtype.type(1.0/2.0))
# 计算 L-4 范数
an = norm(at, 4)
self.check_dtype(at, an)
# 断言 L-4 范数近似为 2^(1/4)
assert_almost_equal(an, an.dtype.type(2.0)**an.dtype.type(1.0/4.0))
# 计算 L-infinity 范数
an = norm(at, np.inf)
self.check_dtype(at, an)
# 断言 L-infinity 范数近似为 1.0
assert_almost_equal(an, 1.0)
# 定义测试向量的方法
def test_vector(self):
# 定义向量 a, b, c
a = [1, 2, 3, 4]
b = [-1, -2, -3, -4]
c = [-1, 2, -3, 4]
# 定义内部函数 _test，用于测试向量的不同范数
def _test(v):
# 断言计算向量 v 的范数近似为 sqrt(30)，小数位数为 self.dec
np.testing.assert_almost_equal(norm(v), 30 ** 0.5,
decimal=self.dec)
# 断言计算向量 v 的 L-infinity 范数近似为 4.0，小数位数为 self.dec
np.testing.assert_almost_equal(norm(v, np.inf), 4.0,
decimal=self.dec)
# 断言计算向量 v 的 L-norm 范数近似为 1.0，小数位数为 self.dec
np.testing.assert_almost_equal(norm(v, -np.inf), 1.0,
decimal=self.dec)
# 断言计算向量 v 的 L-1 范数近似为 10.0，小数位数为 self.dec
np.testing.assert_almost_equal(norm(v, 1), 10.0,
decimal=self.dec)
# 断言计算向量 v 的 L-(-1) 范数近似为 12.0 / 25，小数位数为 self.dec
np.testing.assert_almost_equal(norm(v, -1), 12.0 / 25,
decimal=self.dec)
# 断言计算向量 v 的 L-2 范数近似为 sqrt(30)，小数位数为 self.dec
np.testing.assert_almost_equal(norm(v, 2), 30 ** 0.5,
decimal=self.dec)
# 断言计算向量 v 的 L-(-2) 范数近似为 (205. / 144) ** -0.5，小数位数为 self.dec
np.testing.assert_almost_equal(norm(v, -2), ((205. / 144) ** -0.5),
decimal=self.dec)
# 断言计算向量 v 的 L-0 范数近似为 4，小数位数为 self.dec
np.testing.assert_almost_equal(norm(v, 0), 4,
decimal=self.dec)
# 遍历向量 a, b, c，并分别调用 _test 方法进行范数测试
for v in (a, b, c,):
_test(v)
# 使用指定的数据类型 self.dt 创建数组，并对每个数组调用 _test 方法进行范数测试
for v in (array(a, dtype=self.dt), array(b, dtype=self.dt),
array(c, dtype=self.dt)):
_test(v)
def test_axis(self):
# 测试轴向操作函数
# 比较使用 `axis` 参数与分别计算每行或每列范数的方法。
A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt)
# 对于不同的范数参数进行迭代测试
for order in [None, -1, 0, 1, 2, 3, np.inf, -np.inf]:
# 计算每列的预期范数
expected0 = [norm(A[:, k], ord=order) for k in range(A.shape[1])]
# 断言计算出的每列范数与预期值相近
assert_almost_equal(norm(A, ord=order, axis=0), expected0)
# 计算每行的预期范数
expected1 = [norm(A[k, :], ord=order) for k in range(A.shape[0])]
# 断言计算出的每行范数与预期值相近
assert_almost_equal(norm(A, ord=order, axis=1), expected1)
B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)
nd = B.ndim
# 对于不同的范数参数和轴组合进行迭代测试
for order in [None, -2, 2, -1, 1, np.inf, -np.inf, 'fro']:
for axis in itertools.combinations(range(-nd, nd), 2):
row_axis, col_axis = axis
if row_axis < 0:
row_axis += nd
if col_axis < 0:
col_axis += nd
if row_axis == col_axis:
# 断言当行轴和列轴相同时会抛出 ValueError 异常
assert_raises(ValueError, norm, B, ord=order, axis=axis)
else:
# 计算给定轴的范数
n = norm(B, ord=order, axis=axis)
# 根据轴的组合确定 k_index 的逻辑（仅在 nd = 3 时有效）
# 如果 nd 增加，这部分逻辑需要调整
k_index = nd - (row_axis + col_axis)
if row_axis < col_axis:
# 计算预期的轴范数值
expected = [norm(B[:].take(k, axis=k_index), ord=order)
for k in range(B.shape[k_index])]
else:
expected = [norm(B[:].take(k, axis=k_index).T, ord=order)
for k in range(B.shape[k_index])]
# 断言计算出的范数与预期值相近
assert_almost_equal(n, expected)
# 定义一个测试方法，用于测试保持维度参数的效果
def test_keepdims(self):
# 创建一个形状为 (2, 3, 4) 的 NumPy 数组 A，数据类型由类的成员变量 self.dt 指定
A = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)
# 错误消息格式字符串，用于在 assert_allclose 失败时输出
allclose_err = 'order {0}, axis = {1}'
# 形状不匹配的错误消息格式字符串，用于在 assert_ 失败时输出
shape_err = 'Shape mismatch found {0}, expected {1}, order={2}, axis={3}'
# 检查 order=None, axis=None 的情况
expected = norm(A, ord=None, axis=None)
found = norm(A, ord=None, axis=None, keepdims=True)
# 使用 assert_allclose 检查 found 是否接近于 expected，将其展平以匹配形状
assert_allclose(np.squeeze(found), expected,
err_msg=allclose_err.format(None, None))
# 预期的形状应为 (1, 1, 1)
expected_shape = (1, 1, 1)
# 使用 assert_ 检查 found 的形状是否与 expected_shape 相匹配
assert_(found.shape == expected_shape,
shape_err.format(found.shape, expected_shape, None, None))
# 向量范数。
for order in [None, -1, 0, 1, 2, 3, np.inf, -np.inf]:
for k in range(A.ndim):
expected = norm(A, ord=order, axis=k)
found = norm(A, ord=order, axis=k, keepdims=True)
# 使用 assert_allclose 检查 found 是否接近于 expected，将其展平以匹配形状
assert_allclose(np.squeeze(found), expected,
err_msg=allclose_err.format(order, k))
# 计算预期的形状，将对应轴的长度设置为 1
expected_shape = list(A.shape)
expected_shape[k] = 1
expected_shape = tuple(expected_shape)
# 使用 assert_ 检查 found 的形状是否与 expected_shape 相匹配
assert_(found.shape == expected_shape,
shape_err.format(found.shape, expected_shape, order, k))
# 矩阵范数。
import itertools
for order in [None, -2, 2, -1, 1, np.inf, -np.inf, 'fro', 'nuc']:
for k in itertools.permutations(range(A.ndim), 2):
expected = norm(A, ord=order, axis=k)
found = norm(A, ord=order, axis=k, keepdims=True)
# 使用 assert_allclose 检查 found 是否接近于 expected，将其展平以匹配形状
assert_allclose(np.squeeze(found), expected,
err_msg=allclose_err.format(order, k))
# 计算预期的形状，将对应轴的长度设置为 1
expected_shape = list(A.shape)
expected_shape[k[0]] = 1
expected_shape[k[1]] = 1
expected_shape = tuple(expected_shape)
# 使用 assert_ 检查 found 的形状是否与 expected_shape 相匹配
assert_(found.shape == expected_shape,
shape_err.format(found.shape, expected_shape, order, k))
class _TestNorm2D(_TestNormBase):
# 2维数组的规范化测试类，继承自_TestNormBase
# Define the part for 2d arrays separately, so we can subclass this
# and run the tests using np.matrix in matrixlib.tests.test_matrix_linalg.
# 将2维数组的部分定义为单独的部分，以便我们可以在matrixlib.tests.test_matrix_linalg中使用np.matrix进行子类化和测试运行。
array = np.array # 使用np.array作为数组创建的函数
def test_matrix_empty(self):
# 测试空矩阵的规范化结果是否为0
assert_equal(norm(self.array([[]], dtype=self.dt)), 0.0)
def test_matrix_return_type(self):
a = self.array([[1, 0, 1], [0, 1, 1]])
exact_types = np.typecodes['AllInteger'] # 获取所有精确整数类型的类型码
# float32, complex64, float64, complex128 types are the only types
# allowed by `linalg`, which performs the matrix operations used
# within `norm`.
# `linalg`只允许float32、complex64、float64、complex128类型，这些类型用于执行规范化中使用的矩阵操作。
inexact_types = 'fdFD' # 不精确类型的类型码
all_types = exact_types + inexact_types # 所有可能的数据类型
for each_type in all_types:
at = a.astype(each_type) # 将数组a转换为当前遍历到的数据类型
an = norm(at, -np.inf) # 计算规范化（范数），使用负无穷范数
self.check_dtype(at, an) # 检查计算后的数据类型
assert_almost_equal(an, 2.0) # 断言计算结果接近2.0
with suppress_warnings() as sup:
sup.filter(RuntimeWarning, "divide by zero encountered")
an = norm(at, -1) # 计算规范化，使用-1范数
self.check_dtype(at, an) # 检查计算后的数据类型
assert_almost_equal(an, 1.0) # 断言计算结果接近1.0
an = norm(at, 1) # 计算规范化，使用1范数
self.check_dtype(at, an) # 检查计算后的数据类型
assert_almost_equal(an, 2.0) # 断言计算结果接近2.0
an = norm(at, 2) # 计算规范化，使用2范数
self.check_dtype(at, an) # 检查计算后的数据类型
assert_almost_equal(an, 3.0**(1.0/2.0)) # 断言计算结果接近sqrt(3)
an = norm(at, -2) # 计算规范化，使用-2范数
self.check_dtype(at, an) # 检查计算后的数据类型
assert_almost_equal(an, 1.0) # 断言计算结果接近1.0
an = norm(at, np.inf) # 计算规范化，使用无穷范数
self.check_dtype(at, an) # 检查计算后的数据类型
assert_almost_equal(an, 2.0) # 断言计算结果接近2.0
an = norm(at, 'fro') # 计算规范化，使用Frobenius范数
self.check_dtype(at, an) # 检查计算后的数据类型
assert_almost_equal(an, 2.0) # 断言计算结果接近2.0
an = norm(at, 'nuc') # 计算规范化，使用核范数
self.check_dtype(at, an) # 检查计算后的数据类型
# 需要一个较低的界限以支持低精度浮点数。
# 它们在第7位上有一个偏差。
np.testing.assert_almost_equal(an, 2.7320508075688772, decimal=6) # 断言计算结果接近给定的值
def test_matrix_2x2(self):
A = self.array([[1, 3], [5, 7]], dtype=self.dt)
assert_almost_equal(norm(A), 84 ** 0.5) # 断言计算规范化后的结果接近sqrt(84)
assert_almost_equal(norm(A, 'fro'), 84 ** 0.5) # 断言计算Frobenius范数后的结果接近sqrt(84)
assert_almost_equal(norm(A, 'nuc'), 10.0) # 断言计算核范数后的结果接近10.0
assert_almost_equal(norm(A, inf), 12.0) # 断言计算无穷范数后的结果接近12.0
assert_almost_equal(norm(A, -inf), 4.0) # 断言计算负无穷范数后的结果接近4.0
assert_almost_equal(norm(A, 1), 10.0) # 断言计算1范数后的结果接近10.0
assert_almost_equal(norm(A, -1), 6.0) # 断言计算-1范数后的结果接近6.0
assert_almost_equal(norm(A, 2), 9.1231056256176615) # 断言计算2范数后的结果接近给定的值
assert_almost_equal(norm(A, -2), 0.87689437438234041) # 断言计算-2范数后的结果接近给定的值
assert_raises(ValueError, norm, A, 'nofro') # 断言计算非法Frobenius范数时抛出ValueError异常
assert_raises(ValueError, norm, A, -3) # 断言计算非法范数时抛出ValueError异常
assert_raises(ValueError, norm, A, 0) # 断言计算范数0时抛出ValueError异常
def test_matrix_3x3(self):
# This test has been added because the 2x2 example
# happened to have equal nuclear norm and induced 1-norm.
# The 1/10 scaling factor accommodates the absolute tolerance
# used in assert_almost_equal.
# 创建一个3x3的矩阵A，使用self.array方法，数据类型为self.dt
A = (1 / 10) * \
self.array([[1, 2, 3], [6, 0, 5], [3, 2, 1]], dtype=self.dt)
# 检查A的2范数是否与给定值几乎相等，绝对容差由1/10缩放
assert_almost_equal(norm(A), (1 / 10) * 89 ** 0.5)
# 检查A的Frobenius范数是否与给定值几乎相等，绝对容差由1/10缩放
assert_almost_equal(norm(A, 'fro'), (1 / 10) * 89 ** 0.5)
# 检查A的核范数是否与给定值几乎相等
assert_almost_equal(norm(A, 'nuc'), 1.3366836911774836)
# 检查A的无穷大范数是否与给定值几乎相等
assert_almost_equal(norm(A, inf), 1.1)
# 检查A的负无穷大范数是否与给定值几乎相等
assert_almost_equal(norm(A, -inf), 0.6)
# 检查A的1范数是否与给定值几乎相等
assert_almost_equal(norm(A, 1), 1.0)
# 检查A的负1范数是否与给定值几乎相等
assert_almost_equal(norm(A, -1), 0.4)
# 检查A的2范数是否与给定值几乎相等
assert_almost_equal(norm(A, 2), 0.88722940323461277)
# 检查A的负2范数是否与给定值几乎相等
assert_almost_equal(norm(A, -2), 0.19456584790481812)
def test_bad_args(self):
# Check that bad arguments raise the appropriate exceptions.
# 创建矩阵A和张量B，数据类型为self.dt
A = self.array([[1, 2, 3], [4, 5, 6]], dtype=self.dt)
B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)
# 当ord为'fro'或'nuc'时，使用axis=<integer>或传递1-D数组会导致
# 抛出ValueError异常
assert_raises(ValueError, norm, A, 'fro', 0)
assert_raises(ValueError, norm, A, 'nuc', 0)
assert_raises(ValueError, norm, [3, 4], 'fro', None)
assert_raises(ValueError, norm, [3, 4], 'nuc', None)
assert_raises(ValueError, norm, [3, 4], 'test', None)
# 当ord为除1, 2, -1或-2之外的任何有限数时，计算矩阵范数时，norm应抛出异常
for order in [0, 3]:
assert_raises(ValueError, norm, A, order, None)
assert_raises(ValueError, norm, A, order, (0, 1))
assert_raises(ValueError, norm, B, order, (1, 2))
# 无效的axis值应该抛出AxisError异常
assert_raises(AxisError, norm, B, None, 3)
assert_raises(AxisError, norm, B, None, (2, 3))
assert_raises(ValueError, norm, B, None, (0, 1, 2))
class _TestNorm(_TestNorm2D, _TestNormGeneral):
pass
class TestNorm_NonSystematic:
def test_longdouble_norm(self):
# 非回归测试：longdouble 类型的 p-范数以前会引发 UnboundLocalError。
x = np.arange(10, dtype=np.longdouble)
old_assert_almost_equal(norm(x, ord=3), 12.65, decimal=2)
def test_intmin(self):
# 非回归测试：有符号整数的 p-范数以前的顺序中进行了 float 转换和错误的绝对值。
x = np.array([-2 ** 31], dtype=np.int32)
old_assert_almost_equal(norm(x, ord=3), 2 ** 31, decimal=5)
def test_complex_high_ord(self):
# gh-4156
# 高阶复数测试
d = np.empty((2,), dtype=np.clongdouble)
d[0] = 6 + 7j
d[1] = -6 + 7j
res = 11.615898132184
old_assert_almost_equal(np.linalg.norm(d, ord=3), res, decimal=10)
d = d.astype(np.complex128)
old_assert_almost_equal(np.linalg.norm(d, ord=3), res, decimal=9)
d = d.astype(np.complex64)
old_assert_almost_equal(np.linalg.norm(d, ord=3), res, decimal=5)
# Separate definitions so we can use them for matrix tests.
class _TestNormDoubleBase(_TestNormBase):
dt = np.double
dec = 12
class _TestNormSingleBase(_TestNormBase):
dt = np.float32
dec = 6
class _TestNormInt64Base(_TestNormBase):
dt = np.int64
dec = 12
class TestNormDouble(_TestNorm, _TestNormDoubleBase):
pass
class TestNormSingle(_TestNorm, _TestNormSingleBase):
pass
class TestNormInt64(_TestNorm, _TestNormInt64Base):
pass
class TestMatrixRank:
def test_matrix_rank(self):
# Full rank matrix
# 满秩矩阵
assert_equal(4, matrix_rank(np.eye(4)))
# rank deficient matrix
# 排名不足的矩阵
I = np.eye(4)
I[-1, -1] = 0.
assert_equal(matrix_rank(I), 3)
# All zeros - zero rank
# 全零矩阵 - 零秩
assert_equal(matrix_rank(np.zeros((4, 4))), 0)
# 1 dimension - rank 1 unless all 0
# 1 维度 - 秩为 1，除非全为 0
assert_equal(matrix_rank([1, 0, 0, 0]), 1)
assert_equal(matrix_rank(np.zeros((4,))), 0)
# accepts array-like
# 接受类数组输入
assert_equal(matrix_rank([1]), 1)
# greater than 2 dimensions treated as stacked matrices
# 大于 2 维的数组被视为堆叠矩阵
ms = np.array([I, np.eye(4), np.zeros((4,4))])
assert_equal(matrix_rank(ms), np.array([3, 4, 0]))
# works on scalar
# 对标量也适用
assert_equal(matrix_rank(1), 1)
with assert_raises_regex(
ValueError, "`tol` and `rtol` can\'t be both set."
):
matrix_rank(I, tol=0.01, rtol=0.01)
# 定义一个测试函数，用于测试 matrix_rank 函数在对称矩阵情况下的行列式计算
def test_symmetric_rank(self):
# 断言对单位矩阵求秩应该得到4，hermitian=True 表示输入为对称矩阵
assert_equal(4, matrix_rank(np.eye(4), hermitian=True))
# 断言对全1矩阵求秩应该得到1，hermitian=True 表示输入为对称矩阵
assert_equal(1, matrix_rank(np.ones((4, 4)), hermitian=True))
# 断言对全0矩阵求秩应该得到0，hermitian=True 表示输入为对称矩阵
assert_equal(0, matrix_rank(np.zeros((4, 4)), hermitian=True))
# 创建一个秩亏矩阵 I，将其最后一个元素置为0
I = np.eye(4)
I[-1, -1] = 0.
# 断言对秩亏矩阵求秩应该得到3，hermitian=True 表示输入为对称矩阵
assert_equal(3, matrix_rank(I, hermitian=True))
# 手动设置一个较小的容差值
I[-1, -1] = 1e-8
# 断言对秩亏矩阵求秩应该得到4，hermitian=True，tol=0.99e-8 表示输入为对称矩阵，并且使用指定的容差值
assert_equal(4, matrix_rank(I, hermitian=True, tol=0.99e-8))
# 断言对秩亏矩阵求秩应该得到3，hermitian=True，tol=1.01e-8 表示输入为对称矩阵，并且使用指定的容差值
assert_equal(3, matrix_rank(I, hermitian=True, tol=1.01e-8))
def test_reduced_rank():
# Test matrices with reduced rank
rng = np.random.RandomState(20120714)
for i in range(100):
# Make a rank deficient matrix
X = rng.normal(size=(40, 10))
X[:, 0] = X[:, 1] + X[:, 2]
# Assert that matrix_rank detected deficiency
assert_equal(matrix_rank(X), 9)
X[:, 3] = X[:, 4] + X[:, 5]
assert_equal(matrix_rank(X), 8)
class TestQR:
# Define the array class here, so run this on matrices elsewhere.
array = np.array
def check_qr(self, a):
# This test expects the argument `a` to be an ndarray or
# a subclass of an ndarray of inexact type.
a_type = type(a)
a_dtype = a.dtype
m, n = a.shape
k = min(m, n)
# mode == 'complete'
res = linalg.qr(a, mode='complete')
Q, R = res.Q, res.R
assert_(Q.dtype == a_dtype)
assert_(R.dtype == a_dtype)
assert_(isinstance(Q, a_type))
assert_(isinstance(R, a_type))
assert_(Q.shape == (m, m))
assert_(R.shape == (m, n))
assert_almost_equal(dot(Q, R), a)
assert_almost_equal(dot(Q.T.conj(), Q), np.eye(m))
assert_almost_equal(np.triu(R), R)
# mode == 'reduced'
q1, r1 = linalg.qr(a, mode='reduced')
assert_(q1.dtype == a_dtype)
assert_(r1.dtype == a_dtype)
assert_(isinstance(q1, a_type))
assert_(isinstance(r1, a_type))
assert_(q1.shape == (m, k))
assert_(r1.shape == (k, n))
assert_almost_equal(dot(q1, r1), a)
assert_almost_equal(dot(q1.T.conj(), q1), np.eye(k))
assert_almost_equal(np.triu(r1), r1)
# mode == 'r'
r2 = linalg.qr(a, mode='r')
assert_(r2.dtype == a_dtype)
assert_(isinstance(r2, a_type))
assert_almost_equal(r2, r1)
@pytest.mark.parametrize(["m", "n"], [
(3, 0),
(0, 3),
(0, 0)
])
def test_qr_empty(self, m, n):
k = min(m, n)
a = np.empty((m, n))
self.check_qr(a)
# Verify QR decomposition for empty matrix
h, tau = np.linalg.qr(a, mode='raw')
assert_equal(h.dtype, np.double)
assert_equal(tau.dtype, np.double)
assert_equal(h.shape, (n, m))
assert_equal(tau.shape, (k,))
# 定义一个测试方法，用于测试 QR 分解函数在 mode='raw' 模式下的行为
def test_mode_raw(self):
# 由于因式分解在不同库之间不唯一，结果可能会有所不同，
# 因此无法根据已知值进行检查。功能测试是一种可能性，
# 但需要更多 lapack_lite 中函数的曝光。因此，这个测试的范围非常有限。
# 注意结果是按照 FORTRAN 顺序排列的，因此 h 数组是转置的。
# 创建一个 3x2 的双精度浮点数数组 a
a = self.array([[1, 2], [3, 4], [5, 6]], dtype=np.double)
# 测试双精度浮点数
# 使用 linalg.qr 函数进行 QR 分解，返回 h 和 tau
h, tau = linalg.qr(a, mode='raw')
# 断言 h 的数据类型为双精度浮点数
assert_(h.dtype == np.double)
# 断言 tau 的数据类型为双精度浮点数
assert_(tau.dtype == np.double)
# 断言 h 的形状为 (2, 3)
assert_(h.shape == (2, 3))
# 断言 tau 的形状为 (2,)
assert_(tau.shape == (2,))
# 对 a 的转置进行 QR 分解测试
h, tau = linalg.qr(a.T, mode='raw')
# 断言 h 的数据类型为双精度浮点数
assert_(h.dtype == np.double)
# 断言 tau 的数据类型为双精度浮点数
assert_(tau.dtype == np.double)
# 断言 h 的形状为 (3, 2)
assert_(h.shape == (3, 2))
# 断言 tau 的形状为 (2,)
assert_(tau.shape == (2,))
# 定义一个测试方法，用于测试 QR 分解函数在除经济模式以外的所有模式下的行为
def test_mode_all_but_economic(self):
# 创建两个不同形状的数组 a 和 b
a = self.array([[1, 2], [3, 4]])
b = self.array([[1, 2], [3, 4], [5, 6]])
# 遍历浮点数类型 'f' 和 'd'
for dt in "fd":
# 将 a 和 b 转换为当前浮点数类型 dt 的数组 m1 和 m2
m1 = a.astype(dt)
m2 = b.astype(dt)
# 对 m1 和 m2 分别进行 QR 分解并检查结果
self.check_qr(m1)
self.check_qr(m2)
# 对 m2 的转置进行 QR 分解并检查结果
self.check_qr(m2.T)
# 遍历复数类型 'f' 和 'd'
for dt in "fd":
# 创建复数类型为 1 + 1j 的数组 m1 和 m2
m1 = 1 + 1j * a.astype(dt)
m2 = 1 + 1j * b.astype(dt)
# 对 m1 和 m2 分别进行 QR 分解并检查结果
self.check_qr(m1)
self.check_qr(m2)
# 对 m2 的转置进行 QR 分解并检查结果
self.check_qr(m2.T)
# 定义一个方法，用于检查 QR 分解在不同模式下的行为
def check_qr_stacked(self, a):
# 此测试期望参数 `a` 是一个 ndarray 或其子类，数据类型是不精确的
a_type = type(a)
a_dtype = a.dtype
m, n = a.shape[-2:]
k = min(m, n)
# mode == 'complete' 模式下的 QR 分解
q, r = linalg.qr(a, mode='complete')
# 断言结果的数据类型与输入一致
assert_(q.dtype == a_dtype)
assert_(r.dtype == a_dtype)
# 断言结果的类型与输入一致
assert_(isinstance(q, a_type))
assert_(isinstance(r, a_type))
# 断言 Q 和 R 的形状与预期相符
assert_(q.shape[-2:] == (m, m))
assert_(r.shape[-2:] == (m, n))
# 断言 QR 分解的乘积等于原始矩阵
assert_almost_equal(matmul(q, r), a)
# 创建单位矩阵
I_mat = np.identity(q.shape[-1])
# 扩展单位矩阵以匹配 Q 的形状
stack_I_mat = np.broadcast_to(I_mat,
q.shape[:-2] + (q.shape[-1],)*2)
# 断言 Q 的共轭转置与自身的乘积等于扩展后的单位矩阵
assert_almost_equal(matmul(swapaxes(q, -1, -2).conj(), q), stack_I_mat)
# 断言 R 的上三角部分与 R 相等
assert_almost_equal(np.triu(r[..., :, :]), r)
# mode == 'reduced' 模式下的 QR 分解
q1, r1 = linalg.qr(a, mode='reduced')
assert_(q1.dtype == a_dtype)
assert_(r1.dtype == a_dtype)
assert_(isinstance(q1, a_type))
assert_(isinstance(r1, a_type))
# 断言 Q 和 R 的形状与预期相符
assert_(q1.shape[-2:] == (m, k))
assert_(r1.shape[-2:] == (k, n))
# 断言 QR 分解的乘积等于原始矩阵
assert_almost_equal(matmul(q1, r1), a)
# 创建单位矩阵
I_mat = np.identity(q1.shape[-1])
# 扩展单位矩阵以匹配 Q 的形状
stack_I_mat = np.broadcast_to(I_mat,
q1.shape[:-2] + (q1.shape[-1],)*2)
# 断言 Q 的共轭转置与自身的乘积等于扩展后的单位矩阵
assert_almost_equal(matmul(swapaxes(q1, -1, -2).conj(), q1),
stack_I_mat)
# 断言 R 的上三角部分与 R 相等
assert_almost_equal(np.triu(r1[..., :, :]), r1)
# mode == 'r' 模式下的 QR 分解
r2 = linalg.qr(a, mode='r')
assert_(r2.dtype == a_dtype)
assert_(isinstance(r2, a_type))
# 断言 R2 与 mode=='reduced' 下的 R1 相等
assert_almost_equal(r2, r1)
# 使用 pytest.mark.parametrize 进行参数化测试
@pytest.mark.parametrize("size", [
(3, 4), (4, 3), (4, 4),
(3, 0), (0, 3)])
@pytest.mark.parametrize("outer_size", [
(2, 2), (2,), (2, 3, 4)])
@pytest.mark.parametrize("dt", [
np.single, np.double,
np.csingle, np.cdouble])
# 定义一个测试方法，测试不同尺寸、数据类型的输入
def test_stacked_inputs(self, outer_size, size, dt):
# 使用随机数生成器创建正态分布的 ndarray，并转换为指定数据类型
rng = np.random.default_rng(123)
A = rng.normal(size=outer_size + size).astype(dt)
B = rng.normal(size=outer_size + size).astype(dt)
# 调用 check_qr_stacked 方法，验证 QR 分解的行为
self.check_qr_stacked(A)
self.check_qr_stacked(A + 1.j*B)
class TestCholesky:
@pytest.mark.parametrize(
'shape', [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
)
@pytest.mark.parametrize(
'dtype', (np.float32, np.float64, np.complex64, np.complex128)
)
@pytest.mark.parametrize(
'upper', [False, True])
def test_basic_property(self, shape, dtype, upper):
np.random.seed(1)
a = np.random.randn(*shape)
if np.issubdtype(dtype, np.complexfloating):
a = a + 1j*np.random.randn(*shape)
t = list(range(len(shape)))
t[-2:] = -1, -2
# Compute the product a^H * a or a * a^H based on transpose permutation
a = np.matmul(a.transpose(t).conj(), a)
# Convert 'a' to specified dtype
a = np.asarray(a, dtype=dtype)
# Compute Cholesky decomposition of 'a'
c = np.linalg.cholesky(a, upper=upper)
# Check A = L L^H or A = U^H U depending on 'upper' flag
if upper:
b = np.matmul(c.transpose(t).conj(), c)
else:
b = np.matmul(c, c.transpose(t).conj())
# Set absolute tolerance based on matrix dimensions and dtype precision
with np._no_nep50_warning():
atol = 500 * a.shape[0] * np.finfo(dtype).eps
# Assert the equality of matrices 'b' and 'a' within tolerance 'atol'
assert_allclose(b, a, atol=atol, err_msg=f'{shape} {dtype}\n{a}\n{c}')
# Check if diagonal elements of 'c' (L or U) are real and non-negative
d = np.diagonal(c, axis1=-2, axis2=-1)
assert_(np.all(np.isreal(d)))
assert_(np.all(d >= 0))
def test_0_size(self):
# Define a subclass of np.ndarray with no additional functionality
class ArraySubclass(np.ndarray):
pass
# Create a zero-sized array of integer dtype and cast it to ArraySubclass
a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
# Compute Cholesky decomposition of 'a'
res = linalg.cholesky(a)
# Assert that shapes of 'a' and 'res' are equal
assert_equal(a.shape, res.shape)
# Assert that 'res' has dtype np.float64
assert_(res.dtype.type is np.float64)
# Assert that 'res' is an instance of np.ndarray for documentation purposes
assert_(isinstance(res, np.ndarray))
# Create a zero-sized array of complex dtype and cast it to ArraySubclass
a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
# Compute Cholesky decomposition of 'a'
res = linalg.cholesky(a)
# Assert that shapes of 'a' and 'res' are equal
assert_equal(a.shape, res.shape)
# Assert that 'res' has dtype np.complex64
assert_(res.dtype.type is np.complex64)
# Assert that 'res' is an instance of np.ndarray
assert_(isinstance(res, np.ndarray))
def test_upper_lower_arg(self):
# Explicitly test the 'upper' argument of Cholesky decomposition
a = np.array([[1+0j, 0-2j], [0+2j, 5+0j]])
# Assert that Cholesky decomposition with 'upper=True' is equal to
# the conjugate transpose of Cholesky decomposition with default 'upper=False'
assert_equal(linalg.cholesky(a), linalg.cholesky(a, upper=False))
# Assert that Cholesky decomposition with 'upper=True' is equal to
# the conjugate transpose of default Cholesky decomposition
assert_equal(
linalg.cholesky(a, upper=True),
linalg.cholesky(a).T.conj()
)
class TestOuter:
arr1 = np.arange(3)
arr2 = np.arange(3)
expected = np.array(
[[0, 0, 0],
[0, 1, 2],
[0, 2, 4]]
)
# Assert that the outer product of arr1 and arr2 matches expected
assert_array_equal(np.linalg.outer(arr1, arr2), expected)
# Assert that ValueError with specific message is raised for invalid input arrays
with assert_raises_regex(
ValueError, "Input arrays must be one-dimensional"
):
np.linalg.outer(arr1[:, np.newaxis], arr2)
def test_byteorder_check():
# Check native byte order and assign appropriate character
if sys.byteorder == 'little':
native = '<'
else:
native = '>'
# 遍历数据类型元组，包括 np.float32 和 np.float64
for dtt in (np.float32, np.float64):
# 创建一个 4x4 的单位矩阵，数据类型为 dtt
arr = np.eye(4, dtype=dtt)
# 创建一个新的数组 n_arr，其字节顺序与本地平台一致
n_arr = arr.view(arr.dtype.newbyteorder(native))
# 创建一个新的数组 sw_arr，其字节顺序为小端序并进行字节交换
sw_arr = arr.view(arr.dtype.newbyteorder("S")).byteswap()
# 断言数组 arr 的字节顺序为 '='，即本地平台的默认顺序
assert_equal(arr.dtype.byteorder, '=')
# 遍历线性代数函数元组，包括 linalg.inv, linalg.det, linalg.pinv
for routine in (linalg.inv, linalg.det, linalg.pinv):
# 普通调用，计算 routine(arr) 的结果
res = routine(arr)
# 使用本地字节顺序 n_arr 调用 routine 函数，断言结果与普通调用结果一致
assert_array_equal(res, routine(n_arr))
# 使用字节交换后的 sw_arr 调用 routine 函数，断言结果与普通调用结果一致
assert_array_equal(res, routine(sw_arr))
@pytest.mark.skipif(IS_WASM, reason="fp errors don't work in wasm")
# 标记为跳过测试，如果在 WebAssembly 环境下，原因是浮点错误在 wasm 中不起作用
def test_generalized_raise_multiloop():
# 测试函数：测试通用异常抛出在多重循环中的情况
# 应该会在 ufunc 内部循环的最后一个迭代之前抛出错误
invertible = np.array([[1, 2], [3, 4]])
non_invertible = np.array([[1, 1], [1, 1]])
x = np.zeros([4, 4, 2, 2])[1::2] # 创建一个4x4x2x2的零数组，并取出一部分
x[...] = invertible # 将可逆矩阵赋值给选定的部分
x[0, 0] = non_invertible # 将不可逆矩阵赋值给特定位置
assert_raises(np.linalg.LinAlgError, np.linalg.inv, x)
def test_xerbla_override():
# 测试函数：检查我们的 xerbla 是否成功链接
# 如果没有成功链接，则会调用默认的 xerbla 程序，该程序会向 stdout 打印消息，
# 并且根据 LAPACK 包可能会中止进程
XERBLA_OK = 255
try:
pid = os.fork()
except (OSError, AttributeError):
# fork 失败，或者不在 POSIX 上运行
pytest.skip("Not POSIX or fork failed.")
if pid == 0:
# 子进程；关闭输入输出文件句柄
os.close(1)
os.close(0)
# 避免生成核心文件
import resource
resource.setrlimit(resource.RLIMIT_CORE, (0, 0))
# 下面的调用可能会中止进程
try:
np.linalg.lapack_lite.xerbla()
except ValueError:
pass
except Exception:
os._exit(os.EX_CONFIG)
try:
a = np.array([[1.]])
np.linalg.lapack_lite.dorgqr(
1, 1, 1, a,
0, # <- 无效的值
a, a, 0, 0)
except ValueError as e:
if "DORGQR parameter number 5" in str(e):
# 成功，重用错误码以标记成功，因为 FORTRAN STOP 返回的是成功
os._exit(XERBLA_OK)
# 没有中止，但我们的 xerbla 没有被链接
os._exit(os.EX_CONFIG)
else:
# 父进程
pid, status = os.wait()
if os.WEXITSTATUS(status) != XERBLA_OK:
pytest.skip('Numpy xerbla not linked in.')
@pytest.mark.skipif(IS_WASM, reason="Cannot start subprocess")
@pytest.mark.slow
def test_sdot_bug_8577():
# 测试函数：回归测试，确保加载某些其他库不会导致 float32 线性代数产生错误结果
#
# 在 MacOS 上，存在 gh-8577 可能触发此问题，可能还有其他情况也会出现这个问题
#
# 在单独的进程中进行检查
bad_libs = ['PyQt5.QtWidgets', 'IPython']
template = textwrap.dedent("""
import sys
{before}
try:
import {bad_lib}
except ImportError:
sys.exit(0)
{after}
x = np.ones(2, dtype=np.float32)
sys.exit(0 if np.allclose(x.dot(x), 2.0) else 1)
""")
# 对于每一个不良库 bad_lib，在模板代码中插入 import numpy as np 之前和之后的位置
for bad_lib in bad_libs:
# 构建完整的代码模板，将不良库 bad_lib 插入到 import numpy np 之前
code = template.format(before="import numpy as np", after="",
bad_lib=bad_lib)
# 使用 subprocess 调用 Python 解释器执行该代码
subprocess.check_call([sys.executable, "-c", code])
# 构建完整的代码模板，将不良库 bad_lib 插入到 import numpy np 之后
code = template.format(after="import numpy as np", before="",
bad_lib=bad_lib)
# 使用 subprocess 调用 Python 解释器执行该代码
subprocess.check_call([sys.executable, "-c", code])
class TestMultiDot:
def test_basic_function_with_three_arguments(self):
# multi_dot with three arguments uses a fast hand coded algorithm to
# determine the optimal order. Therefore test it separately.
# 生成随机的 6x2 矩阵 A
A = np.random.random((6, 2))
# 生成随机的 2x6 矩阵 B
B = np.random.random((2, 6))
# 生成随机的 6x2 矩阵 C
C = np.random.random((6, 2))
# 断言 multi_dot([A, B, C]) 的结果与 A.dot(B).dot(C) 几乎相等
assert_almost_equal(multi_dot([A, B, C]), A.dot(B).dot(C))
# 断言 multi_dot([A, B, C]) 的结果与 np.dot(A, np.dot(B, C)) 几乎相等
assert_almost_equal(multi_dot([A, B, C]), np.dot(A, np.dot(B, C)))
def test_basic_function_with_two_arguments(self):
# separate code path with two arguments
# 生成随机的 6x2 矩阵 A
A = np.random.random((6, 2))
# 生成随机的 2x6 矩阵 B
B = np.random.random((2, 6))
# 断言 multi_dot([A, B]) 的结果与 A.dot(B) 几乎相等
assert_almost_equal(multi_dot([A, B]), A.dot(B))
# 断言 multi_dot([A, B]) 的结果与 np.dot(A, B) 几乎相等
assert_almost_equal(multi_dot([A, B]), np.dot(A, B))
def test_basic_function_with_dynamic_programming_optimization(self):
# multi_dot with four or more arguments uses the dynamic programming
# optimization and therefore deserve a separate
# 生成随机的 6x2 矩阵 A
A = np.random.random((6, 2))
# 生成随机的 2x6 矩阵 B
B = np.random.random((2, 6))
# 生成随机的 6x2 矩阵 C
C = np.random.random((6, 2))
# 生成随机的 2x1 矩阵 D
D = np.random.random((2, 1))
# 断言 multi_dot([A, B, C, D]) 的结果与 A.dot(B).dot(C).dot(D) 几乎相等
assert_almost_equal(multi_dot([A, B, C, D]), A.dot(B).dot(C).dot(D))
def test_vector_as_first_argument(self):
# The first argument can be 1-D
# 生成随机的长度为 2 的一维数组 A1d
A1d = np.random.random(2) # 1-D
# 生成随机的 2x6 矩阵 B
B = np.random.random((2, 6))
# 生成随机的 6x2 矩阵 C
C = np.random.random((6, 2))
# 生成随机的 2x2 矩阵 D
D = np.random.random((2, 2))
# 断言 multi_dot([A1d, B, C, D]) 的形状为 (2,)
assert_equal(multi_dot([A1d, B, C, D]).shape, (2,))
def test_vector_as_last_argument(self):
# The last argument can be 1-D
# 生成随机的 6x2 矩阵 A
A = np.random.random((6, 2))
# 生成随机的 2x6 矩阵 B
B = np.random.random((2, 6))
# 生成随机的 6x2 矩阵 C
C = np.random.random((6, 2))
# 生成随机的长度为 2 的一维数组 D1d
D1d = np.random.random(2) # 1-D
# 断言 multi_dot([A, B, C, D1d]) 的形状为 (6,)
assert_equal(multi_dot([A, B, C, D1d]).shape, (6,))
def test_vector_as_first_and_last_argument(self):
# The first and last arguments can be 1-D
# 生成随机的长度为 2 的一维数组 A1d
A1d = np.random.random(2) # 1-D
# 生成随机的 2x6 矩阵 B
B = np.random.random((2, 6))
# 生成随机的 6x2 矩阵 C
C = np.random.random((6, 2))
# 生成随机的长度为 2 的一维数组 D1d
D1d = np.random.random(2) # 1-D
# 断言 multi_dot([A1d, B, C, D1d]) 的形状为 ()
assert_equal(multi_dot([A1d, B, C, D1d]).shape, ())
def test_three_arguments_and_out(self):
# multi_dot with three arguments uses a fast hand coded algorithm to
# determine the optimal order. Therefore test it separately.
# 生成随机的 6x2 矩阵 A
A = np.random.random((6, 2))
# 生成随机的 2x6 矩阵 B
B = np.random.random((2, 6))
# 生成随机的 6x2 矩阵 C
C = np.random.random((6, 2))
# 生成全零的 6x2 矩阵 out
out = np.zeros((6, 2))
# 调用 multi_dot([A, B, C], out=out) 并将结果存储在 ret 中
ret = multi_dot([A, B, C], out=out)
# 断言 out 是 ret 的引用
assert out is ret
# 断言 out 的值几乎等于 A.dot(B).dot(C)
assert_almost_equal(out, A.dot(B).dot(C))
# 断言 out 的值几乎等于 np.dot(A, np.dot(B, C))
assert_almost_equal(out, np.dot(A, np.dot(B, C)))
def test_two_arguments_and_out(self):
# 生成随机的 6x2 和 2x6 的矩阵 A 和 B
A = np.random.random((6, 2))
B = np.random.random((2, 6))
# 创建一个 6x6 的全零矩阵 out
out = np.zeros((6, 6))
# 使用 multi_dot 函数计算 A 和 B 的乘积，并将结果保存到 out 中
ret = multi_dot([A, B], out=out)
# 检查 out 和返回值 ret 是同一个对象
assert out is ret
# 检查 out 的值接近于 A 点乘 B 的结果
assert_almost_equal(out, A.dot(B))
# 检查 out 的值接近于 np.dot(A, B) 的结果
assert_almost_equal(out, np.dot(A, B))
def test_dynamic_programming_optimization_and_out(self):
# multi_dot 函数在传入四个或更多参数时会使用动态规划优化，
# 因此需要单独进行测试
A = np.random.random((6, 2))
B = np.random.random((2, 6))
C = np.random.random((6, 2))
D = np.random.random((2, 1))
# 创建一个 6x1 的全零矩阵 out
out = np.zeros((6, 1))
# 使用 multi_dot 函数计算 A、B、C 和 D 的乘积，并将结果保存到 out 中
ret = multi_dot([A, B, C, D], out=out)
# 检查 out 和返回值 ret 是同一个对象
assert out is ret
# 检查 out 的值接近于 A 点乘 B 点乘 C 点乘 D 的结果
assert_almost_equal(out, A.dot(B).dot(C).dot(D))
def test_dynamic_programming_logic(self):
# 测试动态规划部分的逻辑
# 这个测试直接来自于 Cormen 书中第 376 页。
arrays = [np.random.random((30, 35)),
np.random.random((35, 15)),
np.random.random((15, 5)),
np.random.random((5, 10)),
np.random.random((10, 20)),
np.random.random((20, 25))]
m_expected = np.array([[0., 15750., 7875., 9375., 11875., 15125.],
[0., 0., 2625., 4375., 7125., 10500.],
[0., 0., 0., 750., 2500., 5375.],
[0., 0., 0., 0., 1000., 3500.],
[0., 0., 0., 0., 0., 5000.],
[0., 0., 0., 0., 0., 0.]])
s_expected = np.array([[0, 1, 1, 3, 3, 3],
[0, 0, 2, 3, 3, 3],
[0, 0, 0, 3, 3, 3],
[0, 0, 0, 0, 4, 5],
[0, 0, 0, 0, 0, 5],
[0, 0, 0, 0, 0, 0]], dtype=int)
s_expected -= 1 # Cormen 使用基于 1 的索引，Python 不是。
# 调用 _multi_dot_matrix_chain_order 函数计算 s 和 m
s, m = _multi_dot_matrix_chain_order(arrays, return_costs=True)
# 只有上三角部分（不包括对角线）是感兴趣的部分。
assert_almost_equal(np.triu(s[:-1, 1:]),
np.triu(s_expected[:-1, 1:]))
assert_almost_equal(np.triu(m), np.triu(m_expected))
def test_too_few_input_arrays(self):
# 检查当输入数组过少时是否会引发 ValueError 异常
assert_raises(ValueError, multi_dot, [])
assert_raises(ValueError, multi_dot, [np.random.random((3, 3))])
# 定义一个测试类 TestTensorinv，用于测试 tensorinv 函数的各种情况
class TestTensorinv:
# 使用 pytest 的参数化装饰器，指定多组参数进行测试
@pytest.mark.parametrize("arr, ind", [
(np.ones((4, 6, 8, 2)), 2), # 参数1：4维数组全为1，指数为2
(np.ones((3, 3, 2)), 1), # 参数2：3维数组全为1，指数为1
])
# 测试非方阵处理函数，期望引发 LinAlgError 异常
def test_non_square_handling(self, arr, ind):
with assert_raises(LinAlgError):
linalg.tensorinv(arr, ind=ind)
# 使用 pytest 的参数化装饰器，指定多组参数进行形状测试
@pytest.mark.parametrize("shape, ind", [
# 文档字符串中的示例
((4, 6, 8, 3), 2), # 参数1：形状为 (4, 6, 8, 3)，指数为2
((24, 8, 3), 1), # 参数2：形状为 (24, 8, 3)，指数为1
])
# 测试 tensorinv 函数返回的数组形状是否符合预期
def test_tensorinv_shape(self, shape, ind):
a = np.eye(24)
a.shape = shape
ainv = linalg.tensorinv(a=a, ind=ind)
expected = a.shape[ind:] + a.shape[:ind]
actual = ainv.shape
assert_equal(actual, expected)
# 使用 pytest 的参数化装饰器，指定多组参数进行指数限制测试
@pytest.mark.parametrize("ind", [
0, -2, # 参数1：指数为0和-2
])
# 测试 tensorinv 函数对于超出指数范围的值是否引发 ValueError 异常
def test_tensorinv_ind_limit(self, ind):
a = np.eye(24)
a.shape = (4, 6, 8, 3)
with assert_raises(ValueError):
linalg.tensorinv(a=a, ind=ind)
# 测试 tensorinv 函数的结果是否与文档字符串中的示例一致
def test_tensorinv_result(self):
# 模仿文档字符串中的示例
a = np.eye(24)
a.shape = (24, 8, 3)
ainv = linalg.tensorinv(a, ind=1)
b = np.ones(24)
assert_allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
# 定义一个测试类 TestTensorsolve，用于测试 tensorsolve 函数的各种情况
class TestTensorsolve:
# 使用 pytest 的参数化装饰器，指定多组参数进行测试
@pytest.mark.parametrize("a, axes", [
(np.ones((4, 6, 8, 2)), None), # 参数1：4维数组全为1，axes为None
(np.ones((3, 3, 2)), (0, 2)), # 参数2：3维数组全为1，axes为(0, 2)
])
# 测试非方阵处理函数，期望引发 LinAlgError 异常
def test_non_square_handling(self, a, axes):
with assert_raises(LinAlgError):
b = np.ones(a.shape[:2])
linalg.tensorsolve(a, b, axes=axes)
# 使用 pytest 的参数化装饰器，指定多组形状参数进行结果测试
@pytest.mark.parametrize("shape",
[(2, 3, 6), (3, 4, 4, 3), (0, 3, 3, 0)],
)
# 测试 tensorsolve 函数的结果是否符合预期
def test_tensorsolve_result(self, shape):
a = np.random.randn(*shape)
b = np.ones(a.shape[:2])
x = np.linalg.tensorsolve(a, b)
assert_allclose(np.tensordot(a, x, axes=len(x.shape)), b)
# 定义一个测试函数 test_unsupported_commontype，测试 linalg 是否优雅处理不支持的类型
def test_unsupported_commontype():
# 创建一个浮点16位数组，期望捕获 TypeError 异常，错误信息中包含 "unsupported in linalg"
arr = np.array([[1, -2], [2, 5]], dtype='float16')
with assert_raises_regex(TypeError, "unsupported in linalg"):
linalg.cholesky(arr)
# 定义一个测试函数 test_blas64_dot，测试 64 位 BLAS 下的矩阵乘法
@pytest.mark.skip(reason="Bad memory reports lead to OOM in ci testing")
def test_blas64_dot():
n = 2**32
a = np.zeros([1, n], dtype=np.float32)
b = np.ones([1, 1], dtype=np.float32)
a[0,-1] = 1
c = np.dot(b, a)
assert_equal(c[0,-1], 1)
# 定义一个测试函数 test_blas64_geqrf_lwork_smoketest，烟雾测试 LAPACK geqrf 函数的 lwork 调用，使用 64 位整数
@pytest.mark.xfail(not HAS_LAPACK64,
reason="Numpy not compiled with 64-bit BLAS/LAPACK")
def test_blas64_geqrf_lwork_smoketest():
# 使用 64 位浮点数类型
dtype = np.float64
lapack_routine = np.linalg.lapack_lite.dgeqrf
m = 2**32 + 1
n = 2**32 + 1
lda = m
# 虚拟数组，不需要正确大小，仅供 LAPACK 函数调用
a = np.zeros([1, 1], dtype=dtype)
# 创建一个dtype类型的大小为1的零数组，用于存储工作空间
work = np.zeros([1], dtype=dtype)
# 创建一个dtype类型的大小为1的零数组，用于存储tau参数
tau = np.zeros([1], dtype=dtype)
# 调用lapack_routine函数进行大小查询，获取计算结果
results = lapack_routine(m, n, a, lda, tau, work, -1, 0)
# 确保返回结果中的信息标志为0，表示调用成功
assert_equal(results['info'], 0)
# 确保返回结果中的m字段与输入参数m相等
assert_equal(results['m'], m)
# 确保返回结果中的n字段与输入参数m相等
assert_equal(results['n'], m)
# 将工作空间大小转换为整数，应该是一个合理大小的整数
lwork = int(work.item())
# 确保lwork在2**32和2**42之间，保证大小合理
assert_(2**32 < lwork < 2**42)
def test_diagonal():
# 在这里我们只测试选择的轴是否与数组 API 兼容（最后两个）。`diagonal` 的核心实现在 `test_multiarray.py` 中进行测试。
x = np.arange(60).reshape((3, 4, 5))
actual = np.linalg.diagonal(x)
expected = np.array(
[
[0, 6, 12, 18],
[20, 26, 32, 38],
[40, 46, 52, 58],
]
)
assert_equal(actual, expected)
def test_trace():
# 在这里我们只测试选择的轴是否与数组 API 兼容（最后两个）。`trace` 的核心实现在 `test_multiarray.py` 中进行测试。
x = np.arange(60).reshape((3, 4, 5))
actual = np.linalg.trace(x)
expected = np.array([36, 116, 196])
assert_equal(actual, expected)
def test_cross():
x = np.arange(9).reshape((3, 3))
actual = np.linalg.cross(x, x + 1)
expected = np.array([
[-1, 2, -1],
[-1, 2, -1],
[-1, 2, -1],
])
assert_equal(actual, expected)
# We test that lists are converted to arrays.
u = [1, 2, 3]
v = [4, 5, 6]
actual = np.linalg.cross(u, v)
expected = array([-3, 6, -3])
assert_equal(actual, expected)
with assert_raises_regex(
ValueError,
r"input arrays must be \(arrays of\) 3-dimensional vectors"
):
x_2dim = x[:, 1:]
np.linalg.cross(x_2dim, x_2dim)
def test_tensordot():
# np.linalg.tensordot is just an alias for np.tensordot
x = np.arange(6).reshape((2, 3))
assert np.linalg.tensordot(x, x) == 55
assert np.linalg.tensordot(x, x, axes=[(0, 1), (0, 1)]) == 55
def test_matmul():
# np.linalg.matmul and np.matmul only differs in the number
# of arguments in the signature
x = np.arange(6).reshape((2, 3))
actual = np.linalg.matmul(x, x.T)
expected = np.array([[5, 14], [14, 50])
assert_equal(actual, expected)
def test_matrix_transpose():
x = np.arange(6).reshape((2, 3))
actual = np.linalg.matrix_transpose(x)
expected = x.T
assert_equal(actual, expected)
with assert_raises_regex(
ValueError, "array must be at least 2-dimensional"
):
np.linalg.matrix_transpose(x[:, 0])
def test_matrix_norm():
x = np.arange(9).reshape((3, 3))
actual = np.linalg.matrix_norm(x)
assert_almost_equal(actual, np.float64(14.2828), double_decimal=3)
actual = np.linalg.matrix_norm(x, keepdims=True)
assert_almost_equal(actual, np.array([[14.2828]]), double_decimal=3)
def test_vector_norm():
x = np.arange(9).reshape((3, 3))
actual = np.linalg.vector_norm(x)
assert_almost_equal(actual, np.float64(14.2828), double_decimal=3)
actual = np.linalg.vector_norm(x, axis=0)
assert_almost_equal(
actual, np.array([6.7082, 8.124, 9.6436]), double_decimal=3
)
actual = np.linalg.vector_norm(x, keepdims=True)
expected = np.full((1, 1), 14.2828, dtype='float64')
assert_equal(actual.shape, expected.shape)
# 使用断言来验证两个浮点数的近似相等性，允许误差到小数点后第三位
assert_almost_equal(actual, expected, double_decimal=3)
